(navigation image)
Home American Libraries | Canadian Libraries | Universal Library | Community Texts | Project Gutenberg | Biodiversity Heritage Library | Children's Library | Additional Collections
Search: Advanced Search
Anonymous User (login or join us)
Upload
See other formats

Full text of "Geometric design technique for two-dimensional analog and digital filters."

GEOMETRIC DESIGN TECHNIQUE FOR 
TWO-DIMENSIONAL ANALOG AND DIGITAL FILTERS 



Cosk&n Zumrutkaya 



Thesis 
Z8 



NAVAL POSTGRADUATE SCHOOL 

Monterey, California 






THESIS 




TWO- 


GEOMETRIC DESIGN TECHNIQUE FOR 
DIMENSIONAL ANALOG AND DIGITAL FILTERS 




by 






Coskun Zumrutkaya 






June 19 7 9 




Thesis 


Advisor: S. R. 


Parker 



Approved for public release; distribution unlimited 



UNCLASSIFIED 



SECURITY CLASSIFICATION OF THIS PAGE (Whan Data Entered) 



REPORT DOCUMENTATION PAGE 



I a.£»OAT SUMII* 



READ INSTRUCTIONS 
BEFORE COMPLETING FORM 



2. GOVT ACCESSION NO. 



1. RECIPIENT'S CATALOG NUMBER 



4. T|TlE fend Subtttl*) 



Geometric Design Technique for Two- 
Dimensional Analog and Digital Filters 



5- TYPE OF REPORT a PERIOD COVERED 

Engineer's Thesis; 
June 197 9 



« PERFORMING ORG. REPORT NUMBER 



7. AUTHORS 

Coskun Zumrutkaya 



•. CONTRACT OR GRANT NUMBERS.) 



>. PERFORMING ORGANIZATION NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, California 93940 



10. PROGRAM ELEMENT. PROJECT, TASK 
AREA « WORK UNIT NUMBERS 



II. CONTROLLING OFFICE NAME ANO AOORESS 

Naval Postgraduate School 
Monterey, California 



12. REPORT DATE 

June 19 7 9 



IS. NUMBER OF PAGES 



153 



U. MONITORING AGENCY NAME ft AOORESSf different from Controlling Ottlee) 



IS. SECURITY CLASS, (of thio report) 

Unclassified 



<Sa. OCCLASSIFI CATION/ DOWNGRADING 
SCHEDULE 



16. DISTRIBUTION STATEMENT for thit Report) 



Approved for public release; distribution unlimited. 



17. DISTRIBUTION STATEMENT (o: thm omotrmct entered In Block 30, II different from Report) 



IS. SUPPLEMENTARY NOTES 



IS. KEY WORDS 'Continue on rereree tide it neceeeerr and Identity or block nimtber) 



Analog Filter 
Digital Filter 
Extended Bode approach 



20. ABSTRACT (Continue an referee eide It I 



•ary and Identity ey bleek 



neer) 



A technique for the use of semi-infinite planes to approximate 
the log magnitude versus log frequency characteristics of 
two-dimensional filters is presented (Extended Bode approach) . 
This technique is applied to quarter plane filters and works 
well for separable transfer functions. Other non-separable 
canonic (basic) transfer functions are also studied in terms 
of their planar approximation. The general technique is shown 



DO , r ™ n 1473 

(Page 1) 



EDITION OF 1 MOV «» IS OBSOLETE 
S/N 0103-014-6601 | 



UNCLASSIFIED 



SECURITY CLASSIFICATION OF THIS PAGE fWhan Dote Entered) 



UNCLASSIFIED 



f*eu*wTv claudication q» Tmt mtQtrw*,^ n»i« i.<, K< 



(20. ABSTRACT, Continued) 



to be useful for the insight it provides as well as being a 
simple approach to design. The double bilinear z-transform 
is studied for use with the two-dimensional analog transfer 
functions to convert them to the digital domain. 



DD Form 1473 
„, 1 Jan 73 
S/N 0102-014-6601 



1 Jan 73 UNCLASSIFIED 

N 0102-014-6601 2 iieuAW classification o* ▼»••• natr**™ o«« *«#•*•*) 



Approved for public release; distribution unlimited. 

Geometric Design Technique for 
Two-Dimensional Analog and Digital Filters 

by 
Coskun Zumrutkaya 
Lieutenant /''Turkish Navy 
B.S.E.E., Naval Postgraduate School, 1977 
M.S.E.E., Naval Postgraduate School, 1979 



Submitted in partial fulfillment of the 
requirements for the degree of 



ELECTRICAL ENGINEER 



from the 

Naval Postgraduate School 
June 1979 



ABSTRACT 

A technique for the use of semi-infinite planes to 
approximate the log magnitude versus log frequency charac- 
teristics of two-dimensional filters is presented (Extended 
Bode approach) . This technique is applied to quarter plane 
filters and works well for separable transfer functions. 
Other non-separable canonic (basic) transfer functions are 
also studied in terms of their planar approximation. The 
general technique is shown to be useful for the insight it 
provides as well as being a simple approach to design. The 
double bilinear z-transform is studied for use with the 
two-dimensional analog transfer functions to convert them 
to the digital domain. 



TABLE OF CONTENTS 

I. INTRODUCTION 8 

A. DIGITAL SIGNAL PROCESSING 8 

B. AREA OF INVESTIGATION AND OBJECTIVES 18 

C. PREVIEW OF RESULTS 20 

II. FUNDAMENTALS OF TWO-DIMENSIONAL 

RECURSIVE FILTERING 21 

A. INTRODUCTION 21 

B. DEFINITIONS : 22 

1. Two-Dimensional Signals 22 

2. Two-Dimensional Systems 28 

3. Causality, Separability, Stability 31 

4. Two-Dimensional Difference Equations - 36 

5. Two-Dimensional z-Tranform 39 

6. Two-Dimensional DFT 45 

C. STRUCTURES OF TWO-DIMENSIONAL FILTERS 4 7 

1. One-Dimensional Digital Filter 
Realization 47 

2 . Two-Dimensional Digital Filter 
Realization 53 

D. STABILITY 61 

III. MATHEMATICAL BACKGROUND FROM ANALYTIC 

PLANAR GEOMETRY 69 

IV. CANONIC FACTORS IN CONTINUOUS FREQUENCY 

DOMAIN 87 

A. SEPARABLE FACTORS 87 

B. NON-SEPARABLE FACTORS 10 8 

C. CONCLUSIONS 118 



V. TWO-DIMENSIONAL DIGITAL FILTER DESIGN 
TECHNIQUES 122 

A. INTRODUCTION 122 

B. TWO-DIMENSIONAL DIGITAL FILTER DESIGN 
PROCEDURE BY DOUBLE Z-TRANSFORM 125 

1. Double z-Transform 125 

2. Design Procedure 130 

C. DIGITAL CANONIC FACTORS AND STABILITY 
CONSIDERATIONS 132 

1. Introduction 132 

2. Separable Factors 132 

3. Non-Separable Factors 133 

4. Conclusions 145 

VI. CONCLUSIONS 148 

LIST OF REFERENCES 150 

INITIAL DISTRIBUTION LIST 153 



ACKNOWLEDGMENT 

I am fortunate to have this opportunity to thank the 
many people who have helped in the preparation of this 
thesis. First of all, I am most indebted to my advisor. 
Dr. Sidney R. Parker, for his encouragement, patience, 
guidance and instruction. 

I want to thank the Turkish Navy for the generous 
opportunity to pursue this research. 

I want to thank Mrs. Rosemary Lande for her excellent 
typing. 

Last, but not least, I would like to thank my wife, 
Sumru and my son Devrim, for their patience and understanding 



I. INTRODUCTION 

A. DIGITAL SIGNAL PROCESSING 

Digital signal processing has its roots in 16th century 
mathematics, especially in the fields of astronomy and the 
compilation of mathematical tables. Today it has become 
a powerful tool in a multitude of diverse fields of science 
and technology. The applications of digital signal pro- 
cessing varies from low-frequency spectrum seismology 
through spectral analysis of speech and sonar into the 
video spectrum of radar systems [1] . 

In their book, Theory and Applications of Digital 
Signal Processing [2] , Rabiner and Gold have combined digi- 
tal signal processing theory, with a variety of applications 
ranging from sonar, radar, communication, music, seismic 
and medical signal processing, with digital component 
technology. This technology is the main driving force for 
progress in this field as well as the general area of com- 
puter design. It was also observed by them, that, although 
the formulation of engineering problems is often as vague 
as those of the "softer" sciences (such as anthropology, 
psychology, etc.), the execution of these problems appears 
to depend on greater and greater accuracy and reproducibility 
The capability of digital systems to achieve a guaranteed 
accuracy and essentially perfect reproducibility is very 
appealing to engineers. 



A digital filter is defined in [3] to be a computa- 
tional process in which a sequence of numbers acting as 
input, is transformed into a second sequence of numbers, 
or output digital signal; where the term "digital" implies 
that both time (the independent variable) and amplitude are 
quantized. 

The field of one-dimensional digital filtering, which 
is outlined in figure 1.1, encompasses recursive, non- 
recursive and Fast-Fourier Transform processing. The terms 
recursive and nonrecursive, instead of IIR (Infinite Impulse 
Response) and FIR (Finite Impulse Response) are preferred 
in this thesis. It is shown that recursive processing is 
much more efficient than nonrecursive processing. Stockham's 
method [4] to perform fast convolution, which later became 
known as the FFT method, improved the efficiency of non- 
recursive techniques, so that comparisons in one-dimension 
are no longer strongly biased toward recursive methods [2] . 

One-dimensional recursive digital filter design depends 
strongly on the effective and well developed continuous 
filter design theory. Moreover, the stability analysis of 
higher order one-dimensional recursive realizations can be 
solved by investigation of the root distribution of its 
factored form representation. It is important to realize 
that continuous domain design techniques and factorization 
property (the fundamental theorem of algebra) exists only 
in one-dimension. 



a 



CO Eh 

>H fa 

en Q 



i— i 

Q 



I 



8 



w 

Q 



B 

g 



£ 



us 

I 

H 

fa fa 

o „ 






z fa 



Eh 



CO fa 



fa 



, M 

O fa 

as 



II 



fa 
CD 

9 



fa 

H 
fa 






O fa 






S3 

H Z 



en 



I 

H 

fa 




388 

S l_T 2 

O H H Q 

u S R ^ 

las§ 

O H fa 



< i< H O H O M 
£ .3 2 H fa U Q 



03 

2 
fa 

Eh 
02 

>H 

cn 

fa 

Eh 
H 
U 
H 

Q 

fa 

Z 
O 
H 

Z 
fa 

H 

a 
i 

m 

z 
o 

fa 

o 



fa 

Hi 

> 

fa 
fa 
> 

o 



fa 
fa 

D 



10 



Areas of application of one-dimensional digital signal 
processing are listed in figure 1.1. There are important 
military and specifically naval applications in addition 
to those listed, such as missile and torpedo guidance, 
launcher control system, fire control, combat information 
collection, dissemination and display. 

Digital filtering in several dimensions, which is over- 
viewed in figure 1.2, has gained considerable importance, 
especially for the two-dimensional case. Figure 1.3 presents 
an example of two-dimensional low-pass and high-pass filtering 
Until 196 6, two-dimensional digital filtering was implemented 
by nonrecursive or convolutional techniques. In this method 
the output is the weighted sum of unit sample responses of 
all past input values. The serious disadvantage of the 
convolutional method is the requirements of a very large 
number of arithmetic operations. 

The development of the Fast Fourier Transform (FFT) in 
1966, reduced the number of arithmetic operations consider- 
ably and is used extensively today. Filtering via FFT is 
accomplished by computing the transform of the input func- 
tion, multiplying by the frequency response of the filter, 
and inverse transforming the result. The recursive 
algorithm, which has in general an infinite impulse response 
[3] , constitutes another technique for the realization of 
two-dimensional digital filters. 

Hall [5] compares the amount of computation required for 
the three filtering techniques. He shows that, in general, 
the FFT and recursive algorithms are preferable to the 

11 





02 




T 




*—* 


co 


o 


S 


1-1 


s 


CO 


En 


3 


en 


W 


>H 


2 


co 


H 




Q 


J 


1 


< 


Z 


Eh 




H 


s 


O 


H 


M 




Q 


•r 1 




K 


fa 


O 


O 


fa 




CO 


>H 


2 


05 
O 


£ 


fa 


Eh 


K 


1 


Eh 


CS3 




i 

1 



fa 

OS 
D 
Eh 
U 
D 
« 
Eh 
CO 



CO 
H 
CO 

>• 
< 

z 
< 

fa 

Eh 

< 
Eh 
CO 



05 


a 


ij 


H 


& b 


§£ 


H CO 


a 


z; 


o 


z 




eg 
fa 

Eh 
CO 

>H 

CO 

fa 
< 

Eh 
H 
U 



fa 

<: 
z 
o 

H 

CO 

fa 

S 

H 
Q 
I 

Z 

fa 

o 

fa 

HI 
> 

fa 
> 

o 



CM 



fa 
fa 

CJ 

H 
fa 



12 




Fig. 1.3. a. An example of two-dimensional filtering 
(An original photograph) [2] 



13 




Fig. 1.3.b. 



An example of two-dimensional filtering 
(Low pass filtered version of Fig. 1.3a) 
[2] 



14 



.'.-^■■'/r.-a"' ,\'^y' : ^*-"' '■' ■ \* 



■ * S ' 



Fig. 1 . 3c. 



An example of two-dimensional filtering 
(High pass filtered version of Fig. 1.3a) 
[2] 



15 



nonrecursive one in number of computations and storage 
requirements. Hall demonstrates that the recursive filter 
algorithm constitutes the best method for large data, i.e., 
it is the fastest and cheapest. 

When processing two-dimensional data by recursive filters 
a fundamental problem exists due to inherent feedback, namely, 
the problem of numerical stability. Since, in general, two 
variable polynomials cannot be factored into a product of 
first and second order real coefficient polynomials in each 
of the variables, it is difficult to solve the stability 
problem for two-dimensional recursive filtering. Conse- 
quently, the majority of papers published in two-dimensional 
digital filtering deal with the design of nonrecursive 
filters which are inherently stable [6-10]. 

There are several papers discussing two-dimensional 
recursive digital filters. For example [11] and [12] 
formulate two-dimensional recursive filters by the z- 
transform and linear difference equations, and although 
they investigate problem areas related to stability and 
realization, the majority of problems remain to be solved. 

The most important applications of two-dimensional 
digital filtering of digital data are in picture processing 
and geophysical data analysis. Picture processing can be 
categorized into: 

a) Digital image restoration and enhancement, 

b) Computer pictorial pattern recognition. 



16 



The methods of image restoration and enhancement are 
applied to invert degradations, such as aberrations, 
atmospheric effect, scanning, motion, and to manipulate 
images to improve viewing phenomena experienced by the human 
eye. Image restoration and enhancement has been used with 
great success in biomedicine, i.e., in extraction of quan- 
titative information from x-ray films, chromosome counting, 
measuring the extent of arteriosclerosis from arteriograms, 
in forensic sciences application, i.e., fingerprint image 
enhancement for automatic classification, and in astronomy, 
i.e., removal of turbulence from astronomical photography 
[13]. 

Pictoral pattern recognition by digital computer has 
gained great importance in satellite surveillance and mili- 
tary reconnaissance. The application of two-dimensional 
digital filters permits the separation of different horizontal 
scans on magnetic and gravity maps in geophysics and can 
be used equally well for structural and topographic maps 
or for any other type of factor which is available in the 
format of a planar grid [14] . Although digital filtering 
in several dimensions is gaining importance in medicine, 
i.e., heart volume measurements, in fire control problems, 
and to analyze complex electronic circuitry, there exists 
no general theory concerning structures, analysis and 
design. 

These considerations can be summarized as follows. 
Since one-dimensional recursive digital analysis and design 



17 



techniques depend strongly on the fundamental theorem of 
algebra and on extensive utilization of continuous design 
theory, it represents a special, i.e., a non-generalizable 
field in the area of N-dimensional recursive digital filtering. 
Multi-dimensional digital filtering is well developed in two- 
dimensions but is based, due to the absence of recursive theory 
analysis and design tools, predominantly on nonrecursive 
filtering schemes. Because of their inherent advantages, much 
of the current research is directed towards recursive filtering 
algorithms and techniques. 

B. AREA OF INVESTIGATION AND OBJECTIVES 

One area of investigation of this thesis is to see how the 
semi-infinite straight line approximation technique of one- 
dimensional filters (based upon the log modulus/log frequency 
or decibel vs. log- frequency plots) can be extended to two 
(or higher) dimensional recursive filters. The advantages of 
log-modulus approach is that a transfer function in factored 
form becomes a linear sum of terms when the logarithm of its 
magnitude is taken, and each of these terms can be approxi- 
mated by semi-infinite straight line segments in the one- 
dimensional case, and semi-infinite planes in the multi- 
dimensional case, when expressed in terms of log-frequency. 

The following are specific objectives: 

1. To postulate typical (canonic) factors for 
two-dimensional filters and to investigate their properties 
in the log-modulus/log frequency domain. 



18 



2. To develop an approximation technique for the 
design of two-dimensional recursive filters to meet given 
specifications in the frequency domain. 

The frequency response of two dimensional filters is 
inherently a volume in the space of magnitude, log cj, and 
log CO2 • This volume is the given or desired specification 
of the filter, which in the separable case can be approxi- 
mated by the combination of semi-infinite planes. 

In this thesis we investigate how we can approximate 
this volume by semi-infinite planes as an extension of 
the approximation used in one dimension, i.e., by semi- 
infinite lines. After determining the planar equations 
which perform the given specification of the filter, we 
investigate ways to get a transfer function; first, in the 
continuous frequency domain and second for two-dimensional 
recursive equations in the discrete variable domain. 

The results of this study will classify those charac- 
teristics that can be achieved by cascading simple canonic 
factors, and indicate limitations that are imposed by 
trying to achieve a particular frequency characteristic 
(such as the fan filter) in as simple form as possible. 

The approach taken here is to investigate the problem 
in the continuous domain first, and then study the trans- 
formation from continuous frequency (differential equation) 
domain to discrete time (sampled data) domain. 



19 



C. PREVIEW OF RESULTS 

We have proposed a geometric method to design two- 
dimensional recursive digital filters. This is a new 
method which is relatively simple and easy with respect to 
other design methods. We have investigated properties of 
some typical (canonic) factors and considered their stability 
properties. In Chapter II we go through the basic properties 
of two-dimensional filtering. In Chapter III we provide 
background in terms of planar geometry for the discussion 
in subsequent chapters. In Chapter IV we investigate the 
properties of typical factors and their stability in continu- 
ous variable domain. The new design technique, including 
stability considerations, starting with the continuous 
frequency domain and extending to two-dimensional recursive 
filters in the discrete domain, is presented in Chapter V. 
The results are summarized and some comments for further 
study are made in Chapter VI. It should be noted that the 
results presented apply to quarter plane filters, as defined 
in Chapter V. 



20 



II. FUNDAMENTALS OF TWO-DIMENSIONAL RECURSIVE FILTERING 

A. INTRODUCTION 

There are many signals that are inherently two-dimensional 
for which two-dimensional signal processing techniques are 
required. The most common examples are, of course, images 
in which the two variables are the spatial coordinates. 
Image processing [15], [16] plays an important role in many 
areas of scientific and technical research. Some examples 
are satellite imagery radiography, radar, and biomedical 
images such as acoustical holograms, medical x-rays, and 
electron micrographs [2]. Not all two-dimensional data 
come from an image. In prospecting for gas and oil, a 
common procedure is to monitor the seismic signals produced 
by detonating an explosive charge [17] . An array of geo- 
phones situated along a line radiating in a vertical or 
horizontal direction from the point of explosion is used 
to record these signals. In this case, the two variables 
are distance (from the point of explosion measured along 
the line) and time (the time at which the seismic signal 
reaches a given geophone or, equivalently , a given distance 
from the explosion) . 

Although two-dimensional signals may be processed by 
one-dimensional systems, it is generally preferable to con- 
sider using two-dimensional systems. Many of the basic 
ideas of one-dimensional signal processing may be readily 



21 



extended to the two-dimensional case. There are some 
very important concepts of one-dimensional systems, however, 
that are not directly extentable to two dimensions. It 
is the goal of this chapter to discuss the basic ideas and 
techniques of two-dimensional signal processing and to 
illustrate them into the context of two-dimensional filter 
design. A more detailed discussion of this theory plus 
an overview of recent work in two-dimensional filtering is 
contained in the excellent survey paper by Mersereau and 
Dudgeon [18] , and H. Chang and J.K. Aggarwal [19] . Other 
useful references are the first few chapters of [1] and 
chapter 7 of [2] . 

B. DEFINITIONS 

1. Two-Dimensional Signals 

One-dimensional signals are functions of a single 
variable and two-dimensional signals are functions of two 
integer variables. Consider the two-dimensional sequence 
x(m,n) where m and n are integer variables, shown in figure 
2.1. As in the one-dimensional case, the notation x(m,n) 
is often short hand for a sampled version of a continuous 
two-dimensional signal x(s,t): i.e., 



x(m,n) = x(mT,,nT 2 ) = x(s,t) 



(2.D 

s=mT 1 ,t=nT 2 



One interpretation of x is a function which assigns a 
(generally complex) number to each integer ordered pair 



22 



- 7 y- it/a M 



* 






Fig. 2.1. Pictorial representation for two- 
dimensional discrete array 



X P £$ flit \-$l> J 



23 



(m,n) . The function x is not defined when either or both 

of its arguments (m,n) is not an integer. Such a signal 

will be interchangeably referred to as either a two-dimensional 

array, a two-dimensional sequence or simply as a two-dimensional 

signal . 

Some useful two-dimensional sequences are defined 
below and shown in figures 2.2 to 2.4. These include: 

1. Digital Impulse or Unit Sample 

1 m = n = 
u (m,n) = \ (2.2) 

elsewhere 

2. Digital Step 

( 

I 1 m,n >_ 
u ,(m,n) = \ (2.3) 

elsewhere 

3. Exponential 

m n . n 
a, a~ m,n 

x(m,n) = X (2.4) 

otherwise 

4. Sinusoid (Complex) 

j ((jo^m+oj^n) 
x(m,n) = e -°° < m,n < +=° (2.5) 

As seen above, the two-dimensional step is related to the 
two-dimensional impulse by 



24 



/ 



■+ 



+ 



+-H-4 




-t 



■f- f- -h-h 

■/ -h-h -h +- /-■/- +-/-++- 
-h- -h ■/- f- -H f 4 A -h A 4 +■ 

f +■ 4 + f 4 4 4 -h + + 
4- 4- 4- -hif- 4 4- 4 -h-f- -h 



-4 -h -i-4- f- -f- + -h-h-h 4 i- 
-h -h 4- -h +- -f-f- 4 1-+- j-4- 

+ j-4 4- f -h- 4- ■/- 4 -h -h 
i--l-4--/-4- -/-4-4--h+- 
-h4—h 4-/- 4- 4 +■ 4 4 



I / / > rn. 



Fig. 2.2. Digital Impulse 



25 



- I * * -f 4 

f +■ + -t-_zf- X XXX X 

-h-h + + f XXXXXX 
-h + + 4-4- XX xxxx 
-h-h-h h-h.XX XXXX, 
-h-t- 4- 4- -4 .XXX XXXX 

+.+ + + i .XiXXxxx 

/ / / / / / y ) } y } ) > " 

4--h4-4-4-4-4-4-i--h4- 4- 

-h -h 4- 4- f -4 4 h^4 4 

4--A 4- 4 4 4 4 4 + 4--h 
-h 4 -h-j-h- 4- 4- -h 4- -h 
+ + -h -/ 4 -L- -h 4 4 4 



Fig. 2.3. Digital Step 



26 



-t- -h H- -h -j- -J- -f- ■f- f- -f- 

-*- + + ++• f--*- ff f 
-/- -f- ■+ 4-+ f 4- 4-f- *■ *■ 

-f- + + -h -h f 1 * f- *-* 

■t- -t + -h HfA- J- j- -f- f- + 



m 



4 4 -h -h 4 4 4- + -f- -h -h +■ 
4- + 4 + + 4- 4- + -h +- -h -t- 

+-4-4 -h 4 4- 4 4- 4 4 f- 
-4 4- 7-/ -h -h A -h +~ 4- 



Fig. 2.4. Exponential 



27 



m n 
u_ 1 (m,n) =11 u o (m 1 ,n 1 ) (2.5a) 



m, =-<» m, =-oo 



2 . Two- Dimensional Systems 

A two-dimensional discrete system is defined as 
the unique transformation T that maps an input sequence 
(x(m,n)} and the initial condition sequence {s(m,n)} into 
an output sequence {y(m,n)} 

y(m,n) = T[{x(m,n) },{s(m,n) }] (2.6) 

In signal processing applications, we assume zero initial 
conditions. With this assumption, the above relation may 
be shortened to 

y(m,n) = T[{x(m,n)}] (2.7) 

The system characterized by T is 

(a) linear if 

T[{ax 1 (m,n) + bx 2 (m,n) }] = aT [ {x 1 (m,n) } ] + bT [ {x 2 (m,n) }] 

(2.8) 
for arbitrary constants a and b. 

(b) shift invariant if 

y(m-m 1 ,n-n 1 ) = T [ {x (m-ir^ , n-n 1 ) } ] (2.9) 

for all m, and n, . 

28 



A hardware and software system which is equivalent to the 
transformation T of (2.6) is called a digital filter. The 
class of complex-valued sequences which are square summable 
will play the role of input and output signals of the 
digital filter. Linearity and shif t-invariance are indepen- 
dent properties of a system in that neither property implies 
the other. Thus the system 

T[x(m,n)] = h(m,n) x(m,n) 

is linear but not shift invariant and the system 



2 

T[x(m,n)] = x (m,n) 



is shift invariant but not linear 

For two-dimensional discrete linear systems, the 
impulse response is defined as the response of the system 
at (m,n) to a two-dimensional digital impulse input at 
(m, /n,) with zero initial conditions, i.e., 

h(m,n;m,,n, ) = T [ {u q (m-m 1 ,n-n 1 ) } ] (2.10) 

where u (m,n) denotes the digital impulse. In particular, 
the impulse response of a two-dimensional discrete LSI 
(Linear Shift Invariant) system is 



29 



h (m,n;m, ,n, ) = h (m-m, ,n-n, ) 



or 



h(m,n) = T[{u o (m,n) }] (2.11) 



The response of a two-dimensional LSI* system with zero 
initial conditions is completely characterized by its 
impulse response. In fact using the identity 



x(m,n) = I x(m, ,n , ) u (m-m, ,n-n 1 ) 



l'"l'"o % V "1' 



m 1 ,n 1 =-« 



I x(m-m 1 ,n-n 1 )u (m 1 ,n 1 ) 



m^n^- 00 



we get 



y(m,n) = T[{ x (m, ,n 1 ) u (m-m, ,n-n, ) } ] 



ra r n r 



l x(m 1 ,n 1 )T[{u (m-m 1 ,n-n 1 ) }] 



m 1 /n 1 =- c 



x (m, ,n, ) h (m-m, ,n-n^) 



m l' n l = ~°° ' 



*LSI = Linear Shift Invariance 



30 



y(m,n) = £ x (m-m^n-r^) h (m, ,n,) 



m 1/ n 1 =-< 



= x(m,n) * h(m,n) (2.12) 

where "*" denotes two-dimensional convolution. In other words, 
for LSI systems the basic convolution theorem is valid. 
3. Causality, Separability, Stability 

In discussing various features of two-dimensional 
discrete systems, it is convenient to have spectral nota- 
tion [19] to represent certain regions of the spatial domain. 

1. Let S(a,3) denote a sector with the angular 
interval (a, 3). In polar form S(a,3) may be described by 

S(a,3) = ( (r,6) |r > 0,a < 9 < 3) (2.13) 

where (a, 3) are measured in radians. In particular, we 
give the following notation for the sectors frequently 
encountered: 

S++ = S[0,tt/2]; S-+ = S[Tr/2,ir] 

S— = S[-tt,-tt/2] ; S+- = S[-tt/2,0] (2.14) 

S*+ = S[0,tt) ; S*- = S[-tt,0) 
where "[" and "]" imply including the boundary and ")" implies 



X 



not including the boundary. 



31 



Notice that S++, S-+, S — , and S+- represent first, second, 
third, and fourth quadrants, respectively, and S*+ and S*- 
represent half-planes which are symmetric to each other with 
respect to the origin. (See fig. 2.5.) 

2. For a two-dimensional sequence (x(m,n)}, the 
set 



E v = { (m,n) |x(m,n) ? 0} (2.15) 



is called the support of (x(m,n)}. 

3. Let {h(m,n)} be the impulse response of a two- 
dimensional discrete LSI system, and let E, be the support 
of (h(m,n) }. 

Then, a two-dimensional LSI system is said to be 
causal if E, is the set S++, and semi-causal if E, is the 
set S* + . They are shown in Fig. 2.5. 

A two-dimensional discrete LSI system is said 
to be separable if its impulse response can be factored into 
a product of one-dimensional responses (Fig. 2.6); i.e., 

h(m,n) = h x (m) h 2 (n) (2.16) 

If the equation (2.16) is not satisfied, the filter is said 
to be nonseparable. The advantages of separable filters is 
that two dimensional convolution (2.12) may be carried out 
as a sequence of one-dimensional convolutions. This can be 
seen by rewriting (2.16) as follows. 



32 




fcm 



Fig. 2.5. Support Area of Impulse Response 

(a) Causal 

(b) Semi-causal S. 

(c) Semi-causal S*_ 



33 



X(«r«.) 


H,(-0 




H 2 («.) 











) 



1' 2 



H(z) = H 1 (z 1 )H 2 (z 2 ) 



Fig. 2.6. Separable Digital Systems 



34 



y(m,n) = j» h 1 (m^) h 2 (r^) x (m-m-^n-n-^ (2.17a) 



m, ,n =-<» 



I h 1 (m 1 ) [ I h 2 (n 1 )x(m-m 1/ n-n 1 ) ] (2.17b) 



m =-oo n,=-oo 



I h 1 (m 1 ) a(m-m.,n) (2.17c) 



m l =_c 



where a (m-m, , n) is a sequence of one-dimensional convolutions 
(obtained by evaluating the terms inside the bracket of 
(2.17b) for each fixed value of m, ) . Equation (2.17c) shows 
that y(m,n) may be obtained by a second sequence one-dimensional 
convolution. 

If both the input sequence x(m,n) and the filter 
impulse response h(m,n) are separable, then it is readily 
seen that the output sequence y(m,n) is also separable. 
In this case we obtain the result [using (2.12) and (2.16)] 

00 00 

y(m,n) = [ I h ± Cm^b (m-m 1 ) J I I h 2 (n 1 ) c (n-i^) ] (2.18a) 
m =-« n =-oo 

= a(m)3(n) (2.18b) 



where 



x(m,n) = b(m) c(n) 



35 



V 



A two-dimensional LSI system is said to be (BIBO) 
stable if every bounded input sequence produces a bounded 
output sequence. This condition is satisfied if and only 
if the two-dimensional impulse response of the system is 
absolutely summable, i.e., 

00 

I |h(m,n) | < c° (2.19) 

m, n=-°° 

As in the one-dimensional case, (2.19) can be shown to be a <e- 
necessary and sufficient condition for stability [20]. One 
problem with (2.19) is that it can be quite difficult to 
evaluate for an arbitrary h(m,n). 

4 . Two-dimensional Difference Equations 

As in the one-dimensional case, two-dimensional LSI 
systems can often be described by a two-dimensional linear 
constant-coefficient difference equation relating the output 
y(m,n) of the filter to its input x(m,n). The most general 
form for such a difference equation is shown by 



I a(k 1 ,k 2 )x(m-k 1 ,n-k 2 ) - I b ( ^^ y (m-Jl^n-j^) = 

(k 1 ,k 2 )eR a (i 1 ,l 2 )eR h 

(2.20) 



wh 



ere R and R are finite sets of spatial grid points, 



called the input and output mask, respectively, and {a(k 1 ,k 2 )} 
and (bUwJU)} are the set of constant coefficients that 



36 



characterize the particular filter. Generally, a difference 
equation has a family of solutions as is the case with a 
differential equation. The difference equation is said to 
be recursive if there exists a sequence of computations which 
yields the output sequence serially from the input sequence 
and that portion of the output sequence which is already 
computed/ including the initial condition. An implicit 
ordering of grid points in the spatial domain is assumed for 
the serial computation of the output sequence. 
There are two cases of (2.20): 

1. R , K in the S++ (first quadrant), b(0,0) ^ 
and y(m,n) = for outside of S++ (zero initial conditions), 
then (2.20) represents a causal system. 

2. If R , R. in S*+, b(0,0) ? and y(m,n) = for 

a o 

outside of S*+ (zero initial conditions), then (2.20) 
represents a semi-causal system. 

In both cases we may assume b(0,0) = 1 without loss of 
generality and write (2.20) as 



y(m,n) = £ a (k-^k^ x (m-k-^n-k^ 
(k r k 2 )£R a 



I b(l 1 ,!L 2 )Y(m-l 1 ,n-l 2 ) 

(2.21) 



(l lf l 2 ) eR b -(0,0) 



It is clear that (2.21) is a recursive equation, and the 
implicit ordering, for example may be given as 



37 



{(0,0); (1,0), (0,1), (1,1); (2,0) ,(2,1) , (0,2) ,(1,2); (2,2);...} 

for a causal system, and 

{(0,0), (0,1) , (0,2) ,...; (-1,1) , (0,1) , (1,1) ;...;... 

... (-1,2) , (0,2) , (1,2) ,.. .; . ..} 

for a semi-causal system. 

If {b(£ 1 ,£ 2 )} = {u Q (£ 1 ,£ 2 )} the equation (2.21) 
becomes 

y(m,n) = I a(k 1 ,k 2 )x(m-k 1 ,n-k 2 ) (2.22) 

(k 1 ,k 2 ) £ R a 

so that the output sequence is a weighted moving average of 
the input sequence. The impulse response of (2.22) is 
unique and is obviously {a(k,,k 2 )}. Such filters are called 
finite extent impulse filters (FIR) because the impulse 
response has only finitely many nonzero terms. FIR filters 
obviously are always BIBO stable because a(k,,k 2 ) is always 
absolutely summable. (It has only finitely many nonzero 
terms . ) 

If, on the other hand, the sequence (b(£,,£ 2 )} 
has several nonzero terms, then the impulse response will 
generally have infinitely many terms. We will say that such 
filters are infinite-extent impulse response filters (IIR) . 



38 



Such filters may be BIBO unstable, and so it is necessary 
to devise a stability test for such filters. 

FIR filters have the advantages that they are 
easier to design, that controlling the phase response is 
easier (in particular, it is easy to produce filters which 
have a linear phase response) , and that since the input- 
output relationship is a finite-extent convolution, Fast 
Fourier Transform (FFT) techniques can be used to speed the 
computation of the output sequence. On the other hand, 
recursive (IIR) filters generally require less high speed 
storage and usually require less computation time than FIR 
filters. Comparison of FIR and IIR filters is contained 
in [5]. 

5 . Two-dimensional z-tranform 

A useful mathematical tool for representation of a 
sequence x(m,n) is its two-dimensional z-transform. Given 
a sequence x(m,n), the z-transform of this sequence is 
defined to be: 

< 

CO oo 

Z[{x(m,n)}] 4 X(z 1 ,z 2 ) = I £ x (m,n) z^ (2.23) 

m=-°° n=-°° 

where z, and z 2 are independent complex variables. The most 
useful property of the z-transform is the fact that the z- 
transform of the convolution of two sequences is the product 
of their z-transforms , i.e.: 



39 



Z[{x(m,n) }*{h(m,n) }] = X(z 1 ,z 2 )H(z 1 ,z 2 ) (2.24) 



Therefore, the input-output relationship of an LSI filter 
may be expressed in the z-transform domain as: 

Y(z 1 ,z 2 ) = H(z 1 ,z 2 )X(z 1 ,z 2 ) (2.25) 

where H(z-,,z 2 ), the z-transform of the impulse response 
is called the transfer function. 

A two-dimensional transfer function evaluated on 
the unit circle is called the frequency response of the 
system, i.e., 

jw, jco 9 ju>, jo> 9 jui, joo 9 

Y(e ,e *) = H(e x ,e )X(e x ,e z ) (2.26) 

when a two-dimensional LSI system is describable by the 
difference equation of (2.20), and we take a z-transform >' 
of this, it follows that 



-k -k 2 
[ I a(k 1 ,k 2 )z 1 z 2 ]X(.z 1 ,z 2 ) 



(k x ,k 2 ) eR a 



- [ I b(£ 1 ,£ 2 )z 1 z 2 Z ]Y(z L ,z 2 ) = 

(i lf l 2 ) el^ 



40 



or 



H(z 1 ,z 2 ) = Y(z 1 ,z 2 )/X(z 1 ,z 2 ) 



' k l ~ k 2 
I a(k 1 ,k 2 )z 1 z 2 



(k l' k 2 )eR a 



-A, -l 
I b(£ 1 ,£ 2 )z 1 x z 2 z 

( *■ -i / & ^ ) s^j-. 



= A(z 1 ,z 2 )/B(z 1 ,z 2 ) (2.27) 

where A and B are the z-transf orms of (a(k,,k 2 )} and 
{b (£■•* & 2 ) } respectively. In particular when (2.27) repre- 
sents a causal or semi-causal system, i.e., (2.27) is the 
z-transform of (2.21), we can give explicit descriptions for 

R and R, in each case. For the causal case we define 
a d 

(Fig. 2.5a) 

R a = * ( k i' k 2> I ° i k l i M l' °i k 2 - M 2 } 

(2.28) 

Rj^ = ( (l lf i 2 ) | <_l 1 <_N lf 0£^ 2 - N 2^ 
For the semi-causal case we define (Fig. 2.5b) 

R a = { (k i' k 2 } 1° i k l- r V k 2 = ° ; " M ct i k li M 6'° < k 2 - M 2 } 

(2.29) 

R b = {(£ 1 ,£ 2 )|0<^ ll N 3 ,£ 2 = 0;-N al £ 1 <N a ,Q<£ 2l N 2 } 

41 




HI 

A 

+) 

•H 

+j 

3 

a 
e 
o 
u 

en 

■H 

a; 
o 
c 

<d 

cr 
w 

4-1 
3 

a 

4-> 

a M 

O 0) 

4J 

(U H 

4-1 <44 

5 4J 

C 

144 13 
O <T5 

3 

s cr 

ft 

^4 4J 

n3 ^4 

•H -H 
Q 4-J 



(T3 
CM 



en 

•H 

fa 



42 




C+ 



2 

— 

i 



-4 
1 



jC 
■P 

•H 

T3 
Q) 

-U 

a, 
g 
o 
o 

CO 

•H 

<D 
O 

c 

CD 

a 1 

CO 

-u 
o 
a 

4-> M 



o 



CD 



O 4-1 

<-< 
w 

3 

e to 

(C U 
U I 

03 g 

•H CD 

a co 



IT) 
CN 



CT> 
•H 

fa 



43 



and in both cases b(0,0) = 1 is assumed. The inverse 
z-transform is defined as 



x(m,n) = j j / X(z,,z ) z™ 1 z^ dz,dz„ (2.30) 

(2njr n „ 12 12 12 

C l c 2 

where c^ and c 2 are suitable closed contours in the z, and 
z 2 planes. 

The two-dimensional Fourier transform of a sequence 
is defined as 



DO)., 3u ? 

X(e ,e *) = F[x(m,n)] 



X(z 1 ,z 2 ) 



2 
(z 1 ,z 2 ) eF 



- j ( co-j m+cu^n) 
£ x(m,n)e x * (2.32) 

m r n=-°° 

For any given sequence the set of (z,,z 2 ) for which the 
z-transform converges is called the region of convergence. 
Uniform convergence of the z-transform requires that the 
sequence be absolutely summable, i.e. 

oo 

I |x(m,n) | iz 1 |" m |z 2 |" n < =o (2.32a) 
m, n=-°° 



44 



Similarly, the Fourier transform X(e ,e ) converges 
uniformly to a continuous function of co, and w_ if 



I |x(m,n) | < oo . (2.32) 

m,n=-°° 



A power series of the form of (2.23) is known as a two- 
dimensional Laurent series in the theory of complex functions 
in two variables . 

6 . Two-dimensional DFT 

Let {x(m,n)} be a finite support sequence of the 
size MxN. Then the two-dimensional discrete Fourier 
transform (DFT) of (x(m,n)} is defined as 



M ~ 1N ~ X -^ ( R k i + £ k 2> 

X(k 1# k 2 ) =11 x(m,n)e M L N z (2.33) 

m=0 n=0 



for <_ k. < M-l and <_ k p < N-l. Similarly, the inverse 
discrete Fourier transform (IDFT) is given by 



45 



X(m,n) 



M-l 

- I 



MN 



k,=0 



N-l 
I 

k 2 =o 



x(k lf k 2 )e 



k, k 2 

+:2TT( ir m + ir n) 



(2.34) 



Note that the DFT corresponds to sampling the Fourier trans- 
form of MxN points, i.e., 



X(k x ,k 2 ) 



= x (e 



Ua>i 300. 
/e 



k l k 2 

a) l =27T lT' a) 2 =27T -N- 



(2.35) 



]W. 



30). 



More generally, if X(e ,e ) is the Fourier transform 
of (x(m,n)} whose size is greater than MxN, then the 
IDFT of {X(k 1 ,k )} will result in aliasing. In other words, 
the IDFT of (X(k,,k 2 )} will yield one period of (x(m,n)} 
where 



x(m,n) = 



m 1 ,n 1 =-<» 



x (m+m,M,n+n,N) 



(2.36) 



46 



C. STRUCTURES OF TWO-DIMENSIONAL FILTERS 

1. One-dimensional Digital Filter Realization 

One-dimensional digital filters, in general, can 
be shown in the z-transform domain as follows: 

y 



1=0 Y(z i ) 

H < z l> = ~ = xuJt (2 ' 37) 

1 - I a m Z l 
m=l 



or , in the time domain 



M L 

y(n) I a y(n-m) + J b x(n-l) (2.38) 



m 
m=l 1=0 



The realization implementing equation (2.38) directly, is 
known as the direct form 1 and shown in Figure 2.7. 

Equation (2.37) can be rewritten in a slightly 
different form by introducing a new variable, W(z 1 ) , such 
that 



W(z,) Y(z,) 

H / z ) = L_ . . x s (2.39) 

H(Z 1 J X(z x ) W z 1 



The corresponding set of difference equations consists of: 



47 



I 



N 



J t— 1 



2 t 3 

I 
H 



as 




kl 



M 




r 



e 

o 

t^ 

+j 

o 
cu 

u 

•H 
Q 



CN 



0^ 
•H 
fa 



o 



5< 



48 



M 
w(n) = x(.n) + I a m w(n-m) (2.40) 

m=l 



and 

L 
y(n) I b £ w(n-£) (2.41) 

1=0 

Equations (2.40) and (2.41) are realized in the 
direct form 2 as shown in Figure 2.8 which can be visualized 
as a cascade realization of two digital filters, realizing 
the denominator and numerator polynomials, respectively. 
Figure 2.8 can be redrawn as shown in Figure 2.9. The 
resulting realization is called direct form 3 or canonic 
form, since it represents a structure with the number of 
delays equal to the order. 

At this point, an additional realization will be 
introduced, the direct form 4, which combines the charac- 
teristics of the direct forms 1 and 2 having only two 
entries in each summer, which corresponds directly to 
hardware implementation and by realizing the numerator and 
denominator polynomials separately in cascaded form 
(Figure 2. 10) . 

There are several other methods to realize a one- 
dimensional M order digital filter. For example, the 
cascade form of first and second order sections H(z,), as 
shown in Figure 2.11, where 



49 



H(z X ) 



2 b £ Z l 
£=0 



-£ 



M 



1 " I a z. 

L ml 

m=l 



-m 



X(n) 




W(n) 



I 



«— 2 


H 


«, r 


H 


rf 1 
^ 3 V 


7-1 




Fig. 2.8. Direct Form 2 



50 



H(z 1 ) 



£=0 * L 



M 



i - y a z" 

L ml 
m=0 



-m 




Y(-) 



Fig. 2.9. Direct Form 3 
(Canonic) 



51 



I 



X t~o II 



J ^ 



N 




WH 



l^H- 




u-jK 





JO 




cv^~~ 


■ ^ 


a 


CjW- 


3 

— 4 








_ — 

-± 

KT 



c 
X 









£ 

O 
fa 

-P 

o 

CD 
S-i 
•H 
Q 



o 

H 

•H 
fa 



52 



H(z 1 ) = a Q n H i (z 1 ), * = ^ (2.42) 
i=l 



Each component H^(z^) can be realized in one of the above 
outlined forms. 

Partial fraction expansion methods applied to 
equation (2.37) lead to a parallel form realization of 
first and second order sections (Figure 2.12), i.e., 

k 
H(z 1 ) - c+ J H i (z 1 ), k= [^+1] (2.43) 

i=l 

Other forms of realization include hybrid structures, 
i.e., parallel-cascade forms, the transpose configuration, 
which can be obtained for all previously outlined structures 
by reversing the direction of signal flow and by interchanging 
all branch and summing modes, wave-digital and continued 
fraction expansion filters. 

The direct forms discussed above and the series/parallel 
arrangement of lower order sections are by far the most 
often used ones in computer simulation and in digital 
hardware. 

2. Two-dimensional Digital Filter Realization 

Two-dimensional Quarter plane digital filters, in 
general form, can be shown as 



53 




N 




ii 

n" 

x 



s 
u 



fa 

0) 

o 
en 

u 



CM 



•H 
fa 



54 




HW = C+2H } U,) 



Fig. 2.12. Parallel Form 



55 



£ =0 £ =0 X Z 
H(z,,z ) = — - (2.44) 

M l M 2 

-m.. -m 

1 ~ I I a mm z, z 
L *< m.iru 1 2 

m l =0 m 2 = ° 

The numerator of the two-dimensional transfer function (2.44) 
can be written as a one-variable polynomial in z~, weighted 
by coefficients which are functions in z,, i.e., 



L 2 


L l 


I 


( I 


«,,=0 


J..=0 


L 2 






' S V 



£ 1 _il 2 



z 2 (2.45) 



" £ 2 
:„ (2.46) 



The same reasoning applied to the denominator polynomials 
leads to 



M M 

-m, -itu 

D ( z r z 2 } = 1 - ^ ( l ^2 Zl )Z2 

m,=0 m.=0 

(2.47) 

M 2 

= 1-1 (V. < z l»» Z 2 m (2 - 48) 
m 2 =0 



56 



To obtain the direct form 3 representation, a new function, 



W(z 1 z 2 ) is introduced: 



W(z,,z 9 ) Y(z,,z 9 ) 

H(z i z 2> - xi^Tzf) ■ w(^fr (2 - 49) 



where each factor can be written using equation (2.46) as 



W(z x ,z 2 ) = 



*2 



£ 2 -0 



/ N I 1 £ 2 (Z 1 ))Z 2 *<V*2 J 

(2.50) 



and equation (2.48) as 



M 2 



-m 



2 



Y(z 1 ,z 2 ) = W(z 1 ,z 2 ) + I (D M m (z 1 ))z 2 Z Y(z 1 z 2 ) 

m 2 =0 ^ 2 

(2.51) 

The realization of equations (2.50) and (2.51) is shown in 
Figure 2.13 using compact notation, and Figure 2.14 where 

N L l ( z l^ anc ^ D M m ^ z l^ are ^ m Pl emented f° r each (L-,^ 2 ) 
and (nunu) by one-dimensional direct form 4 realizations. 

The general form of two-dimensional causal digital 
filter (2.44) can be rewritten by introducing a new variable, 
W(z-,z 2 ) or in the previous function, then, 

W ( z , z « ) Y ( z , , z 9 ) 

H <V Z 2> " X( Zl ,z 2 ) • W(z',z 2 ) (2 " 52 ' 



57 




e 
o 

•H 

-p 
+J 



2 

-u 
o 

(d 

a 
e 
o 
u 



a 



e 
o 

+J 

u 
aj 
u 

■H 

Q 



m 

H 

•H 



58 




Q 



CO 

e 

S-l 


Pn 

-P 
U 
<U 

u 

•H 

Q 



■H 

fa 



59 



Each factor can be written as : 



L l L 2 

-I -I 

Y(l l'»2 J = I I W Z l 1 z 2 2 WCz r z 2 ) 

£ 1 =0 £ 2 =0 



(2.53) 



and 



M M 

-m -m« 
W(z r z 2 ) =- XC V z 2 ) + I I a z x z 2 W(z r z 2 ) 

m 1 =0 m 2 =0 - 1 - 



(m 1+ m 2 )^0) (2<54) 



The equations (2.53) and (2.54) can be written as a differ- 
ence equation as follows, 



L l L 2 
y(m,n) = I I b z z W(m-2. 1 ,n-£ 2 ) (2.55) 

l l =0 £ 2 =0 1 2 



and 



M l M 2 

w(m,n) = x(m,n) + I I a m m w ^ m " m 2 ' n ~ m 2 * (2.56) 

m 1 =0 m 2 =0 



m 1 ,m 2 ^ 



60 



Without loss of generality, L. is chosen equal to M. for 

all i. The structure realizing equation (2.55) and (2.56) 

is shown in Figure 2.15. It is noted that the number of 

delays in Figure 2.15 equals the order of H(z..,z ). The 

structure is, therefore, by definition, canonic. 

The direct form 4 (Fig. 2.10) realization of H(z.,z 2 ) 

is shown in Figures 2.16 and 2.17, using compact notation and 

one-dimensional direct form 4 (Fig. 2.10) implementation of 

each N « and D , respectively. 
L 1^2 M l m 2 

D. STABILITY 

Two major problems in the design of a recursive filter 
are approximation and stability. Since the output values 
are used through feedback by the recursive algorithm, it is 
possible for the output values to become arbitrarily large 
independent of the input. A filter of this kind is said to 
be unstable, which is an undesirable condition. Thus we 
need to know what constraints to put on the recursive filter 
coefficients, so 'that the filtering operation will be stable. 
This problem is similar to the stability problem for the one- 
dimensional case, except that the added dimension inherently 
increases the complexity of the analysis. The major diffi- 
culty is that the concept of zeros and poles does not hold 
in multi-dimensional systems so that simple algebraic systems 
extrapolation of one-dimensional results is not obvious. 
• Stability Theorems : 

Theorem 1 (Shank's theorem): A causal recursive filter 
with the z-transform H(ZwZ 2 ) = A(z 1 , z 2 ) /B (z 1 , z 2 ) , where A 



61 




(Q*— I — hQ 



m 






(hh— I — ►© 



M 



Qk 



4 



© 



t-j 













TCNJ 




c 



c 

U 



Q 

CN 



CO 

e 
o 

-U 

u 

CD 
U 
•H 
Q 



U") 



CN 



•H 

fa 



62 




rc\j 

Tsl 



og 



— @* — 









M 





rsj + 






i 
x 



c 
o 

■H 
+J 

id 

-u 



c 

+J 
o 
td 

cu 

2 

o 



Q 



u 
o 

4J 

o 

0) 
S-l 
•H 
Q 



^o 



CN 



Cn 
■H 



63 



rsj 






H 

% 

J5' 




Q 
CM 



6 
S-i 

o 
fa 

u 

0) 



CM 



t7> 
■iH 

fa 



and B are polynomials in z, and z 2 , is stable if and only 
if there are no values of z, and z such that B(z ,z„) = 
for all z. _> 1 and z« >^ 1. (In general we assume 

B(Z 1' Z 2 } = b 00 + ^O 2 !" 1 + b 01 Z 2" 1 + ^l 2 !* 1 ^" 1 ••" 
that is an expression in increasing powers of z. and z 9 .) 

In other words, the theorem says that if there are any 
values (real or complex) of z, and z 2 for which B(z,,z ? ) 
is zero for z. and z 2 simultaneously greater than or equal 
to unity in magnitude, then the filter will be unstable. 

It is much more difficult to determine the stability of 
two-dimensional filters than of one-dimensional filters. 
In the one-dimensional case, it is only necessary to locate 
finite set of roots in the z-plane. In order to mechanize 
the theorem, we must, in general, find the values of z. and 
z_ for which B(z,,z~) = 0. There is no technique for locating 
the zeros of a general two-dimensional polynomial in (z,,z_). 
Shanks, Tretial and Justice, in their paper [11] proposed a 
technique to apply theorem 1 to determine the stability of 
a two-dimensional recursive filter, H(z,,z 2 ), by which the 
unit disk in the z, plane must be mapped into the z^ plane 
by solving the implicit two-variable denominator polynomial 
of H(z,,z 2 ). A necessary and sufficient condition for 
stability is determined if the map of the unit disk 
d, = (z-,z_ _> 1) does not overlap the unit disk 
d 2 E (z ? ,z > 1) in the z plane. This method requires an 
infinite number of mappings and, therefore, can not be 
applied exactly. 



65 



Huang [20] simplified Shank's stability theorem 
considerably by stating the stability theorem as follows 
Theorem 2: A causal filter with a z-transform 



A(z, ,z.J 
H(z w z n ) = 



'l'"2' B(z 1 ,z 2 ) 



where A and B are two-variable polynomials in z,,z , is 
stable if and only if: 

1) the map of the unit circle 6d = (z , | z, | =1 ) , 
according to B(z ,z 2 ) = 0, lies outside of the unit disk 
d_ = (z ,|z |_>1) in the z -plane, and 

2) no point in the unit disk d.. = ( z ^ , | z , | _> 1 ) 
maps into the point z = by the relation B(z..,z 2 ) = 0. 

In order to apply this theorem as a stability test, we 
have to map the unit disk 5d, = (z , |z 1 |=1) into the z ? 
plane according to B(z, ,z 2 ) = or z 2 = f (z,) | i i-, and 
to see whether the contour in the z^ plane lies inside the 
unit disk d 2 = (z- , | z 2 |<_1) . Also, it is necessary to solve 
B(z,,0) = to see whether there are any roots with magnitude 
greater than 1. 

Although theorem 2 is much simpler than using Shank's 
original theorem, the procedure is still infinite in the 
sense that, in principle, we have to map the unit circle 
5d, into the z ? -plane. 

Ansell [21] , whose main contribution is to couple the use 
of a Hermite test for checking stability with a series of 



66 



Sturm tests checking positivity, has reduced theorem 2. 
Although Ans ell's results enable us to test stabilty in a 
finite number of steps, it still is, unfortunately, very 
tedious. 

Make the change of variables: 






and, 



1 - Z 2 

1 + z 2 



and let 



E(p, /P?) 
H(ZwZ ) = 



l'~2 



fcp 1 ,p 2 ; 



where E and F are polynomials in p, and p~. We can restate 
Theorem 2 as follows. 

Theorem 3 (Ansell's theorem): The causal recursive 
filter H(z,,z 2 ) is stable if and only if: 

1) for all real finite w , the complex polynomial 
in p„, F(joj,p 2 ) has no zeros in ReP 2 > 0, and 

2) the real polynomial in p-,, F(p,,l) has no 
zeros in Rep, 0. 

This theorem is essentially the same as Theorem 2, but 
an advantage of Theorem 3 is that condition 1) can be 
tested using standard techniques of circuit theory. 



67 



Anderson and Jury [22] modified Huang's method by out- 
lining a procedure which replaces the bilinear transforma- 
tions by the construction of a Schur-Cohn matrix and checking 
for positivity of a set of self-inversive polynomials. It 
also replaces the Hermite test by a Schur-Cohn matrix test 
and requires a series of Sturm tests or equivalently , a 
system of tests establishing the roots distribution of a 
polynomial. These modifications represent a substantial 
reduction in computations as compared to Huang's method. We 
will discuss stability considerations which are related to 
our design technique in Chapter IV and Chapter V in detail. 

In this chapter, basic properties of two-dimensional 
digital filters have been reviewed and summarized including 
semi-causality to establish a fundamental theory of two- 
dimensional digital filters. 



68 



III. MATHEMATICAL BACKGROUND FROM ANALYTIC PLANAR GEOMETRY 

In this study, we propose using semi-infinite planes to 
approximate a given frequency response of a two-dimensional 
analog filters by extending one-dimensional semi-infinite 
straight line approximation to two-dimensional cases. For 
example, as we will see, the frequency response (dB vs. log- 
frequency) of the separable canonic factor 

H(ja) 1 ,j(u 2 ) = cl+jUiTi ) (l+ja) 2 T 2 ) 

can be approximated by a collection of planes. This is 
discussed in Chapter IV in detail. 

To prepare for the discussions in the next chapters, 
it is worthwhile to review and develop useful equations 
from analytic planar geometry which are going to be used 
throughout this study with some examples of how to use them. 

The properties of the analytic planar geometry can be 
summarized as follows: 

1. The general equations for a plane are given by 

Ax+By+Cz+D=0 (3.1) 

where A, B and C are the direction numbers of the normal 
to the plane. The equation of a plane passing through the 
point (xwy^z,) and perpendicular to the line which has 
direction numbers of A, B and C is given by: 



69 



A(x-x 1 ) + B(y-y 1 ) + CU-z^ = (3.2) 

2. The direction cosines of the normal to the plane 
are given by- 



cos a = — ^ = T-TTo (3.3a) 

(A Z +B Z +C ) 1/Z 



cos 6 = — ? ^ — TT7T (3.3b) 

(A z +B +C ) ±/Z 



cos y = 



{A 2 + B 2 + C 2 ) l/2 



where a, 8/ and y are the angles that the normal makes with 
the x, y, z axes respectively. 

The relationship between direction cosines is given by 



2 2 

Cos a + Cos 3 + Cos y = 1 (3.3d) 



3. For a line joining points (x, ,y^/Z-^) and (x 2 ,y2> z 2) 
the equation is given by 



Cos a _ Cos 3 _ Cos y ,^ ^, 

x 2 -x x y 2 -y! z 2 -z 1 



4. The equation of a straight line with the direction 
numbers a, b, c and passing through the point (x^y^z^ is 
given by 



x- Xl y-y 1 z- Zl 

a ~ b = n — (3.5) 



5. The angle 6, between two intersecting lines with 
direction numbers a, b, c and a', b f , c' respectively is 
given by: 

Cos = Cos a Cos a 1 + Cos 6 Cos 3' + Cos y Cos y 1 

(3.6a) 
or 



Cos 9 = aa ' + bb ' + cc ' - (3.6b) 



Va 2 +b 2 + c 2 Va ,2 +b ,2 + c'2 



6. For parallel lines, the relationship between 
direction numbers is given by 



a' b^ c* U * /J 



7. For perpendicular lines: 

aa' + bb' + cc ' = (3.8) 

8. The direction numbers of the line of the intersection 
between two given intersecting planes, with direction numbers 
A, B, C and A' , B', C respectively, is given by 



71 



a = 



b = 



B 



B* C 



A 



C A' 



(3.9a) 



(3.9b) 



c = 



B 



A' B' 



(3.9c) 



The angle 6 between these two planes is given by 



Cos 



AA' + BB 1 + CC* 



V 



A 2 +B' 



+ C V A ' 



(3.10) 



2 2 2 



9. The equation for a plane in terms of its slope can 
be calculated as follows: 

As seen from Fig. 3.1, the slope, m, of the normal 
line is given by 



m = tan 



h 

g 



^hFW 



(3.11) 



m = slope of the plane 



M = tan (90° - 9) 



= - 2. = . 1 
h m 



(3.12) 



72 



normal to olane 




normal to nlane 



plane 



Fig. 3.1. A General Plane Configuration 



73 



The equation of the plane is given by one of the 
following: 

A(x-x 1 ) +By + Cz = (3.13a) 

Ax + B(y-y,) + Cz = (3.13b) 



and 



Ax + By + C(z-z 1 ) = (3.13c) 



A*l = By ± = Cz 1 = -D (3.13d) 



The intersection of the plane with the z-x plane 
is given by 



A(x -x 1 ) + Cz = (3.14) 



and the slope, nu, of this line is 



»i = a§ = -§ (3 - 14a) 



Similarly, the intersection of the plane with the 
z-y plane is given by 

B(y -y, ) + Cz = (3.15) 



74 



and the slope of this line, nu , is 



m 2 ~ " dy~ " " C (3.15a, 



Substituting (3.15a) and (3.14a) into (3.13) 
yields the following alternate equations for the plane. 

-m(x-x-j) - m 2 y + z = (3.16a) 

-m,x - m 2 (y-y 2 ) + z = (3.16b) 

-itux - m 2 y + Cz - Zj) = (3.16c) 

The direction cosines of the normal in terms of 
m-, / m 2 are 

-m, 

Cos a = (3.17a) 



Vm^ 



2 2 



-m~ 
Cos 6 = (3.17b) 



v^7 



m 2 2 + l 



Cos Y = 



Vm-L 2 + m 2 2 +1 



(3.17c) 



Substituting (3.14a) and (3.15a) into (3.11) and 
(3.12) yields 



m = i- ( 3 - 18a) 



m 2 



75 



and 



M = -V m i +m 2 (3.18b) 



2 



The angle, 8, between two intersecting planes from 
(3.10), (3.14a) and (3.15a), is given by 



m,m, * + irunu ' + 1 ,- ,„, 

Cos 9 = LJ: Li (3.19) 



ymj 2 + m 2 2 + 1 Vmp+mp + 1 



The direction numbers of the line of intersection, 

from (3.9), (3.14a) and (3.15a) are given by 

a = (m£ - m 2 ) CC (3.20a) 

b = (m£ - m 1 ) CC* (3.20b) 

c = (m|m 2 - m^ir^) CC * (3.20c) 

The direction cosines of this intersection line 
are given by 



m 2 ~ m 2 (3.21a) 

Cos a = 



V2 2 

(ml - nu) + (m| - m, ) + (m-jnu - mlm^) 



2 



m' - m, 
Cos 6 = (3.21b) 



V2 2 

(mi -nu) + ^ m i ~ m i) + (m^m 2 -m^ir^) 



76 



m l m 2 " m 2 m l ,, ,. , 

Cos y = (3.21c) 



V^ 7 



2 2 2 

m 2 ) + (m* -m,) + (m|m - m'm, ) 



10. When two planes are added in one dimension to form 
a new plane, the following results. Given two planes 



-m,x - nuy + z - z, = 



-m,'x - mly + z 1 - z' = 



Their sum in the z-direction with z" = z* + z is 
given by 

-(m 1 +m')x - (m 2 +m 2 ')y + z" - ^ + z£) = 

The slopes of the intercepts with the zx and zy planes 



m 1 " = -(m 1 + m|) 



m 2 " = -(m 2 + m 2 ') 



The slope of the normal is given by 

2 , ,_ , .,2,-1/2 



m 



" = [ (m, +m[) + (m 2 +m 2 ) ] 



(3.22 



77 



To aid in comprehension, two examples follow: 
Example (1) : Given m, (slope in z,x plane), m, 
(slope in x,y plane) , and x, (the point at which x-axis goes 
through the plane). Find the equation of the plane. 
The general plane equation (3.1) is given by 

Ax + By + Cz + D = 

When z = 

B D 
x A Y " A 



m 3 = If = " I (given) 



and 



When y = 



and 



x, = - j (given) 



A D 
" " C x " C 



■l - E - "I (given) 



D 
Z l C 



78 



If we substitute these into the general equation, the result is 



or 



A 

Ax - Am-,y - — z - Ax, = 
3 J m, 1 



x - nuy - — z - x, = 
3 m, 1 

or 



z = nux - nunuy - x i m i (3.23) 

As shown here, we can derive the plane equation easily 
from the given specification. After finding this equation, 
we can find the canonic factor by simply letting 



x = log go. 



y = log a)- 

i i 2 
z = 10 log |T| 

Alternately, we can go from the canonic factor to equations. 
Thus equation (3.23) can be written as 

2 

10 log |T(joo 1 , joo 2 ) | = m 1 log oj 1 - m^ log oj 2 - 10 log K 

(3.24) 
where 

10 log K = ^i x i 



79 



This result assumes that 

i i2 w l 

|T(ju>,,ju> 2 ) | = (3.24a; 

Taking the logarithm of (3.24a) yields 

2 

10 log |T| = 10 p log u^ - 10 q log oj 2 - 10 log K 



(3.25) 



By equating (3.25) and (3.24), we get 



10 p = m 1 (3.26a) 



10 q = m 1 m 3 (3.26b] 



It follows from (3.26a) and (3.26b) that 



m 3 = q/p . 



Thus the transfer function of (3.24a) has a logarithmic 
characteristic as sketched in Fig. 3.2. 

The previous result can be extended to the following 
transfer function 



(1 + j Wl ) p/2 

T(j<D, / ja> 9 ) = —t-t? 7TT? (3.27) 

1 Z K ±/Z (l + jw 2 ) q/ 



By substituting (3.26a) and (3.26c) into (3.27), then 



80 



Z=io 1^1T| 




■s loovJ 



*=: looi*/| 



Fig. 3.2. The logarithmic characteristic of 

P 



T(jgo 1 , jo) 2 ) 



jj- 



Kcu. 



20 log |t| = m, log co, - m, m log aj_ - 10 log K 



81 



nu/20 

(1 + jco,) L 

T( >l' jw 2 ) ■ — Hu/20 (3 * 28) 

K 1/Z (l + ja) 2 ) Z 



When w, << 1 and oj 2 << 1 , region I 
T(jo Jl /jco 2 ) = K~ 1/2 

and 

20 log |T| = - 10 log K 

When co, << 1 and oj- >> 1* region II 



T(j Ul ,ja> 2 ) = — ^-720 

K 1/2 (ja J2 ) A 



and 



20 log |T| = - 10 log K - nu log cj. 



When oj 2 << 1 and oj, >> 1, region III 

111,720 

(j«2.) 

T(ju) 1 ,jo) 2 ) = -j/2 

K 



and 



20 log |T| = - 10 log K + m, log 



0). 



82 



When u>, >> 1 and oo- >> 1, region IV 



nu/20 

j OJ- 



T(ju), ,jw~) = , 

1 2. . nu/20 



K X / 2 H 1 2 ' 



and 



20 log |T| = - 10 logK + m 1 log oj, - m, log to,. 

The resulting planar approximation for the logarithmic 
transfer function is shown in Figure 3.3 
Example 2 ; 

In this example, we derive planar equations from the 
given transfer function. 

Given a separable two-dimensional transfer function: 

T(jw, , ju 9 ) = ■= (3.29) 

(1 + JU) 1 T 1 ) P (1+ JU) 2 T 2 ) H 

The magnitude square of the transfer function (3.29) is 
given by 

K 2 

|T(ja) 1 ,ja) 2 )| = 2 2~p" 2 2~q < 3 - 3 °) 

1 A (1 + o Ji Z t 1 Z ) P (1 + ^ 2 T 2 } 

Take the logarithm of both sides of (3.30), 



10 log |T| 2 = -10 p log(l+co 1 2 x 1 2 ) - 10 q log(l+u> 2 x 2 ) + 20 log K 



(3.31) 



when w,!, << 1 and co 2 x 2 << 1, the equation (3.31) becomes 



83 



A fcrtoUjITl 




t 7^ta>W 




IE 



Br 



3- l°j^ 



> TUlaj^l 



Fig. 3.3. The logarithmic characteristic of 

(1 + jo)) Mi/20 

X Z K ±/Z (l + jco 2 ) il 2/ zu 

20 log |t I = -10 log K + m, log u>, - m- log go- 
84 



2 

10 log |T| = 20 log K = C 



where C is some constant. 

When w^t^ >> 1 and oj 2 t 2 >> 1, 



or 



10 log |T| = -20p log u^ - 20q log a, T 



+ 20 log K 



i i 2 
10 log |T| = -20p log o^ - 20q log a> 2 - 20p log t. 



Let 



-20q log x 2 + 20 log 



K 



2 

z = 10 log |t| 



x = log to. 



y = log a). 



and we get the planar equation, 



Z = - 



20px - 20qy - 20p log t, - 20q log t 2 + 20 log K 



or 



85 



p q 

T T 

(-20p)x + (-20q)y - z = 20 log (-± ? ) 



K 



or 



n^x + m 2 y - z - v 



where 



m^ = slope of the line (x,z) plane 
= -20 p db/decade. 

m 2 = slope of the line (y,z) plane 
= -2 0q db/decade. 

In a similar manner the regions (co, t, << 1; 

CJ 2 T 2 >> "^ an< ^ ^l 1 ! >> ^' (J °2 T 2 <<c ^ can ^ e i nvest: '- < ? ate( i' 
The planar approximation for the logarithmic transfer function 

is sketched in Fig. 4.1 where q = p = 1. 

In this chapter we have developed some useful planar 

geometry formulas and shown how they can be used to determine 

a planar approximation for the two-dimensional logarithmic 

transfer function for separable cases. Some non-separable 

factors are discussed in the next chapter in which experimental 

verification of this approximation technique is also presented. 



86 



IV. CANONIC FACTORS IN CONTINUOUS FREQUENCY DOMAIN 

In Chapter III, we discussed some equations of planar geo- 
metry and showed how to use them for approximating separable two- 
dimensional transfer functions. In order to improve the scope 
of the proposed design technique, we now consider the behavior 
of other separable and non-separable factors insofar as 
their planar geometry approximation is concerned. This will 
provide insight into the frequency response characteristics 
of two-dimensional systems and indicate what can be expected 
by cascading these canonic factors, knowing their individual 
behavior. The logarithmic characteristic of cascaded factor 
is simply the sum of their individual logarithmic character- 
istics. This approach parallels the semi-infinite line 
approximation (dB vs. log w ) technique which has been used 
so successfully in one-dimensional design. 

A. SEPARABLE FACTORS 

The following several separable canonic factors are now 
considered. 

(1) (1 + jaj 1 T 1 ) P (l + jt 2 u3 2 ) q (4. a) 

i 

(2) (1 + jaJ 1 T 1 ) P (l + ju) 2 T 2 ) q (l + joo 2 T 3 ) q (4.b) 



(3) [1 + (ja Jl T 1 ) P ] a • (4.c) 



where p, q and a are integers. 

87 



(1) First, as an example of the technique, consider 
the following analog separable transfer function, the first 
separable canonic factor with p = q = -1; 

H(ju> ,ju> ) = ,-, 7 ■ w-i ■ • r (4.1) 

1 2 (1 + joj 1 t 1 ) (1 + joo 2 t 2 ) 

The logarithm of magnitude-square of this function 
is given by 

*> *j _ 

10 log |H(co 1 ,oo 2 ) | = -10 logd+u^ t^ 2 ) - 10 log(l+co 2 t 2 ) 

(4.2) 

Let us define the following for simplicity of notation. 

z = 10 log |H(oo 1 ,o) 2 ) I = 20 log |H(o) 1 ,a) 2 ) | (dB) 

(4.3a) 
x = 20 log oo-l (4.3b) 



y = 20 log w 2 



(4.3c) 



Consider the regions in the (gj, ,u) 2 ) plane (as seen 
in Figure 4.2) as follows: 

Region I, when oo,x, < 1 and ^ 2 x 2 < 1, 

then z = db (4.4a) 

Region II, when co^t -^ > 1 and co 2 t 2 < 1 

z = -20 log (w 1 t 1 ) db (4.4b) 



88 



Region III, when ^^l < 1 and W 2 T 2 > 1 

z = -20 log(co 2 T 2 ) db (4.4c) 



and 



Region IV, when ou^ > 1 and u t~ > 1 



z = 



= -20 log (oj 1 t 1 ) - 20 log (oj 2 t 2 ) (4.4d) 

The above equations (4.4a), (4.4b), (4.4c) and (4.4d) 
can be expressed as 

Region I, z = (4.5a) 

Region II, z = -x - 20 log t-, (4.5b) 

Region III, z = -y - 20 log t 2 (4.5c) 

Region IV, z = -x - y - 20 log t 2 - 20 log x, 

(4.5d) 

These equations describe four intersecting planes in 
x,y,z space as shown in Figure 4.1. 

Plane I is the horizontal plane at zero db. Plane II 
and Plane III have slopes of -20 db/decade passing through 
the corner lines as indicated in Figure 4.1. These lines 
are analogous to the break point frequencies of one-dimensional 



89 



Z = 20 I o g |H( w r w 2 ) I 




og w 



slope —-20 d B 



- 20dB 



Fig. 4.1. Planar Approximation for 

1 



H (co, ,u>~) 



1+J0J 1 T 1 ) (1 + J00 2 T 2 ) 



2 2 



20 log I H (oi 1 ,co 2 ) | = -10 log (1+w, T, 

-10 log (1+co 2 2 t 2 2 ) 



90 



*2* 1 /^ 



Y = I O 



g w 



IT 




W 



I 






X = 



°g w. 



Fig. 4.2. The regions of 



H(j<VJ0) 2 ) " (1 + J03 1 T 1 ) (1 + J0) 2 T 2 ) 



2 2 



20 log |H(o) 1 ,o) 2 ) | = -10 log (1+u^ Tj ) 

2 
2 T 2 



2 2 

-10 log (1+oj„ t 9 ) 



in (log u), , log co 9 ) plane 



91 



log-modulus plots. Several properties of this plane are 
summarized. Pertinent equations from analytic geometry are 
given in Chapter III. 

The actual frequency response of this separable 
canonic factor with t, = t, = 1 is shown in Figure 4.3a 
(dB vs. log w,, log oj 2 ) which fits the planar approximation , 
In Fig. 4.3b/ the contour plot of the frequency response is 
given, which verifies the regions in log uj-, / log u>- plane 
shown in Fig. 4.2. The logarithmic plots versus linear 
frequency are given in Fig. 4.3c and Fig. 4.3d. It is 
apparent that the approximation regions are well defined 
when logarithmic frequency is used. 

The frequency response of this separable canonic 
factor with t, = 1, t- = 2 is given in Figures 4.4. These 
results again confirm the logarithmic approximation. 

(2) Next consider the canonic form 

H(j Ml ,ja> 2 ) = (1 + j^tjj (1 + ju) 2 x 2 ) (1 + ja> 2 T 3 ) 

where we have two separable factors of oj 2 . The frequency 
response of this factor with t, * 1, x 2 = 2 and x^ = 4 is 
shown in Figures 4.5. As seen from this figure, we have 
-40 dB/decade slope in the direction of u 2 , as expected 
from the one-dimensional case. The regions are shown in 
the contour plot with dotted line (Fig. 4.5c) . 



92 



In a similar manner we can analyze these regions as 
follows: The logarithm of magnitude square of this function 
is given by 



10 log |H(u) 1 ,u) 2 ) | 2 = - 10 log(l + oj 1 2 t 1 2 ) - 10 log (1 + w 2 2 x 2 2 ) 



- 10 log(l + w 2 2 t 2 ) 



Consider the following regions in (to, , w 2 ) plane (as 
seen in Fig. 4.5a) as follows: 

Region I, when u-,t-, << 1 and co 2 x 2 < ^ 2 t 3 << 1 
then z = dB. 

Region II, when u,t, >> 1 and o) 2 t 2 < ^t-d <k 1 

z = - 20 log (oj- l t 1 ) dB 
Region II, when co-,t-, << 1 and oj 2 t 2 < ^2 T 3 <<; 1 

z = -20 log (co 2 x 2 ) - 20 log(oj 2 x 3 ) dB 
Region IV, when m,T, >> 1, ^2 T 3 > C0 2 T 2 >:> ^ 



z = -20 log u,t, - 20 log ^ 2 T 2 " 20 l0g W 2 T 3 



93 



Region V, when co^ << 1, ^^^ >> 1, ^ T >>o) 2 t 2 >> 1 



Z = -20 logU 2 T 3 ) - 20 log(o, 2 T 2 ) 



Region VI, when ^^l >> ^' W 2 T 3 > T 2 (J °2 >:> 1 



Z = - 20 log a) 1 T t - 20 log oo 2 t 2 - 20 log oj 3 t 3 



where 



Z = 20 log |H(u) 1 ,t 2 ) | . 
(3) Consider the separable canonic factor now, 



H(a) 1 ,o) 2 ) = [1 + (jw 1 T 1 ) P ] a 



We can show that the following regions occur by 
using the notation of (4.3a) and (4.3b) for simplicity: 



94 





Fig. 4.3a. The frequency response of 

. 1 

H(jo) 1 ,:o) 2 ) - (i +jWi ) (l+jo) 2 ) 

2 

z = 2Q log|H(o) 1/ o) 2 ) | = -10 logd+u^ ) 

(dB vs. log bu,log o> 2 ) 



-10 log (l+w 2 ) 



95 




X-SCflLE=l .OOE+00 UNITS INCH 
T-SCRLE=I . OOE+00 UNITS INCH 



•ig. 4.3b. The contour plot of H^,^) = (1+j ) ( i +j(0 ) 

z = 20 log |H| = -10 log (l+o^ 2 ) - 10 log (l+io 2 2 ) 



96 





Fig. 4.3c. 



The frequency response of 
H(joJ 1 ,jco 2 ) 



(l+jai 1 ) (l+jo) 2 ) 



z = 20 log|H| = -10 log(l+oj 1 2 ) -10 logd+o^ 2 ) 
(dB vs. 03, ,0)-) (linear frequency) 



97 



10. 



>", 



a-, 



3. 



I 

T 




Fig. 4.3d. The contour plot of H(oo..,go_) = . . , . , . = , . r 

r 1 2 (1+]U) 1 ) (l+:co 2 ) 

z = 20 log|H| = -10 logd+u^ 2 ) - 10 log(l+w 2 2 ) 

(linear frequency) 



98 




Fig. 4.4a. The frequency response of 
H(do) 1 ,do) 2 ) 



(l+joa 1 ) (l+j2co 2 ) 



2 2 

= 20 log|H| = -10 log(l+oj 1 ) - 10 log(l+4co 2 ) 

(dB vs. log co, , log co 2 ) 



99 




X-SCRLE-1 .OOE+00 UNITS INCH. 
T-SCflLE=l . OOE+00 UNITS INCH. 



Fig. 4.4b. 



The contour plot of the frequency response of 
H(o3 1 w 2 ) = 



(l+ja) 1 ) (l+j2u) 2 ) 



z = 20 log|H| = -10 log(l+co 1 2 ) - 10 log(l+4to 2 2 ) 



in (log to, , log to 2 ) plane 



100 








Fig. 4.4c. The frequency response of 

1 



H ( j co , , j co -, ) = ,-, , • - 

J 1 2' (l+jco,) (1 + j2co 2 ) 

z = 10 log|H| = -10 logd+co, 2 ) - 10 log(l+4co 2 2 
(dB vs. co, ,co_) (linear frequency ccale) 



101 




Fig. 4.4d. The contour plot of the frequency response of 



H(jo) 1 ,joj 2 ) = 



(1+w, ) (l+4oo 2 ) 
z = 20 log|H| = -10 log(l+oo 1 2 ) - 10 log(l+4aj 2 ) 
(dB vs. o) 1 ,o) 2 ) (linear frequency scale) 



102 





^ ' o g w 2 




«H 


"2 


YE 




nr 


12 


0)i*Tx 








i 


I 






UH-i 



-► !og w 1 



Fig. 4.5a. 



The region of 
H(j(D 1 ,ju) 2 ) = 



(l+jo) 1 T 1 ) (1+joo t ) (l+jco 2 x 3 ) 
in (log a), /log co„) plane 



103 




Fig. 4.5b. The frequency response of 

1 
H(jo) 1 ,jo) 2 ) - (i+j 0Jl ) (l+j2co 2 ) (l+34w 2 ; 

z = 20 log|H| = -10 log(l+03 1 2 ) - 10 log(l+4w 2 ) 

2 
-10 log(l+16u> 2 ) 

(dB vs. log o) lf log w 2 ) 



104 




X-SCRLE=1 . OOE+00 UNITS INCH. 
Y-SCRLE=I . OOE+00 UNITS INCH. 



Fig. 4.5c. The contour plot of the frequency response of 
H(ju lf ju 2 ) = ( i+j Ul ) (l+j2oo 2 ) (l+34co 2 ; 

2 2 

z = 20 log|H| = -10 log(l+a) 1 ) - 10 log (1+4^ ) 

2 
-10 log(l+16oj 2 ) 

(log oo-L vs. log co 2 ) 



105 




Fig. 4.5d. 



z = 



The frequency response of 

H(joJ 1 ,jOi 2 ) = (1+j(Ui) (i +j 2co 2 ) (l+j4oo 2 ) 

20 log|H| = -10 log(l+ca 1 2 ) - 10 log(l + 2aj 2 2 ) 

-10 log(l+16oo 2 2 ) 
(dB vs. (ji-.,uiy) (linear frequency scale) 



106 



10, 



0. - 



X- 



/«u 




Fig. 4.5e. The contour plot of the frequency response of 

. 1 

H(D(o 1 ,30) 2 J - (1+jc0i ) (l+j2oa 2 ) (l + j4co 2 ) 

? 2 

z = 20 log|H| = -10 log (1+w^) - 10 log(l + 4aj 1 ) 

-10 log (l+16oo 2 2 ) 
(linear frequency scale) 



107 



Region I, when go, > 1/t, , 

z = apx + 20ap log t-, dB 
Region II, when to, < 1/t-, 

z = dB 

These separable canonic factors are identical to 
the canonic factors of a one-dimensional transfer function. 
They are just an extension to the two-dimensional case. 
Therefore, they need not be discussed further. 

B. NON-SEPARABLE FACTORS 

The following non-separable canonic factors are now 
considered. 

(1) [1+ (jaj 1 x 1 ) P (joj 2 T 2 ) q ] a (4.d) 

(2) [1 + (ju 1 T 1 + joo 2 T 2 ) P ] a (4.e) 

where p, q and a are integers. 

(1) First, consider the non-separable canonic 
factor of the form: 



H(oJ 1 ,aj 2 ) = [1 + (jaj- L T 1 ) P (joJ 2 T 2 ) ] 



108 



Region I, when 



P q 1 

1 2 T P T q 



z = dB 



Region II, when 



P q 



1 2 T p , q 



T l T 2 



z = apx + aqy + 20ap log t, + 20 ay log t- dB 



The cornerline between the two regions is given by the 
equation 



or 



p q p q 

u l w 2 T l T 2 



px + qy + 20p log t, + 20q log t 2 = 



The cornerline is sketched in Figure 4.6. 

The slope of the plane in region II is given by 



-4 2 



'p + q . Direction numbers of the normal are given 
by -ap, -aq and -1. 



109 



»,= « 2 = I 



(*)(*!* 



^ x - lo g 



I 



T 



M±)H^ 



y -lo g 



Fig. 4.6. The regions of non-separable case in 
log to, , log co^ plane 



110 



Depending upon the sign of p and q and the magnitude of 
t, and x 2 several different orientations of the corner line 
are possible. 

Example (1) 

Discuss this canonic form in the special case with 
p = q = 1, a=+l, x, =0.5, and x 2 = 1. The transfer 
function gets the following form: 

H(o),,a) 2 ) = [1 + (.5jw,) ( ja) 2 ) ] 

Take the logarithm of the magnitude square. 

i 2 | i+2 

z = +20 log|l - .5(jj,co 2 | 
In region I, w,t, <<: ^ an< ^ aj ? T 2 <<: "*" 



z = dB. 



In region II, id,t, >> 1 and w 2 t 2 >> 1 
z=x + y+20 log (0.5) 



where 



x = 20 log co 1 , y = 20 log cu 2 and z = 20 log |H(oj 1 ,u 2 ) 

The actual frequency response of this factor is shown in 
Fig. 4.7 including a contour plot. 

Ill 




Fig. 4.7a. 



The frequency response of 

H( joo 1 , joo 2 ) = 1 + ( ja>,) ( jO. 5o3 2 ) 

z = 20 log|H| = 20 log (1 - 0.5oj,oj 2 ) 

(dB vs. log go,, log oo^) 



112 




X-SCRLE=l;OOE+00 UNITS INCH. 
T-SCRLE=1..0£iE + 00 units inch. 

Fig. 4.7b. The contour plot of the frequency response of 
H(jco L ,ja> 2 ) = 1 + ( jco 1 ) (j0.5oj 2 ) 
z = 20 log|H| = 20 log (1 - O.Su^u^) 
in (log co, , log ai 2 ) plane 



113 



When a = -1 we have a stability problem in this example 
The denominator of this transfer function becomes zero 
when 



or 



1 — . 5o)-j 0)2 = 



o)-| cop — 2 



This is the singularity of the transfer function in (w,,^) 
plane- It can be called a singular line of the system. It 
is a curve in the (ui,,to 2 ) plane to be compared with a pole 
in a one-dimensional transfer function. 

Example #2 ; 

In this example, we discuss the same canonic form 
with different coefficients, i.e., 

p = l, q = -1 , a = +l, 

x, = 0.5, x 2 = 1 

Then the transfer gets the form of 

* 5uj l +1 
H( Ul ,« 2 ) = (1 + ~— > 



114 



Taking the logarithm of its magnitude squared 
|H(a), ,u 9 ) = ll + ±\ Z 

. 5(jd, 

z = +20 log ll + 



w 2 



. 5a).. 

In Region I when - << 1 

a) 2 



z = dB. 

• 5ai, 

In Region II, when >> 1 

w 9 



z = 20x - 20y + 20 log (.5) 

where x = log u^, y = log u 2 and z = 20 log |h(u>,/u> 2 ) |. The 
boundary between these regions is given by the equation 

. 5oi-| = oi~ 

We also have a stability problem for this canonic factor when 
u>2 s 0. This is also a singularity in (u), ,w 2 ) plane in the 
form of a straight line. This type of factor is unstable. 
(2) Now, consider the non-separable canonic factor 



H(o),,a> 2 ) = [1 + (ju),w, + ju^oO ] 



We have two regions again, as follows: 
Region I, when w,t, + w 2 t 2 < 1 

z = dB 

115 



Region II, when u),x, + w-x- > 1 

z = ap log (oa 1 x 1 + aj 2 T 2 ) • 

The separation line between the two regions is 
given by the equation: 



U 1 T 1 + 0) 2 T 2 = 1 



or 



log (10, x, + o) 2 t 2 ) = . (4.6) 

The locus of points of (4.6) is shown in Figure 4.8, 
For region II, we have to think of two conditions, 
as follows: 

a) When w,x, > oo 2 x 2 

z = apx + ap20 log x, 

b) When oj,t, < ^ 2 x 2 

z = aqy + ap20 log x 2 . 

Thus, region II is split into two areas with the separation 
line given by 



W 1 T 1 = W 2 T 2 ' 



116 




Fig. 4.8. Locus points of logC^T^ + w 2 x 2 ) - 



117 



Frequency, response of this non-separable transfer function 
is shown in Fig. 4.9 with p - 1, a = -1, x = x = l. 

C. CONCLUSIONS 

In this chapter, we discussed the properties of the 
separable and non-separable two-dimensional properties of 
canonic factors in the analog domain (co,,^ plane). We can 
conclude that these canonic factors can be approximated by 
planes in the log magnitude versus log co, and log oj^ space. 
As a design technique we can obtain any given frequency 
specification by proper combinations of planar approximation 
of canonic factors. It is obvious that only certain special 
cases and symmetries can be obtained using the factors dis- 
cussed in this thesis. For example, a wedge-shaped character- 
istic in frequency domain could not be synthesized by means 
of the factors we have considered. 

The other result of this study is the stability problem 
which arises in certain transfer functions, that is, the 
singularity curve or line. For a separable case, we don't 
have a stability problem; but in non-separable cases we may 
have stability problems. 



118 





Fig. 4.9a. 



The frequency response of 
J 1 J 2 1 + joo, + jog. 



z = 20 log|H| = 10 log 1 1+ (oj, +oj 2 ) 
(dB vs. log to, /log co 2 ) 



-, 2 



119 




X-SCRLE=1 . OOE + 00 UNITS INCH 
T-SCflLE-1 . OOE+00 UNITS INCH 



Fig. 4.9b. 



The contour plot of the frequency response of 
H(ju>, ,jw ) = . —-. 

z = 20 log|H| = 10 log |l + (u^+u^) 2 ! 
in (log co, , log co-) plane 



120 





Fig. 4.9c. The frequency response of 

1 
H(3oa 1 ,:a) 2 ) = 



1 + ju). 



JO). 



z = 20 log|H| = 10 log |l + (u^+a^) | 
(dB vs. a) 1 ,oo 2 ) (linear frequency scale) 



121 






v - TWO-DIMENSIONAL DIGITAL FILTER DESIGN TECHNIQUE 

A. INTRODUCTION 

The design of two and multi-dimensional recursive digital 
filters is difficult due to the absence of the Fundamental 
Theorem of Algebra in two or more dimensions [15]. That is, 
polynomials in two variables may not, in general, be factored 
into lower order polynomials. For this reason, several design 
procedures for two-dimensional recursive filters have been 
proposed in the literature. We then consider the use of the 
double bilinear z-transfer with the semi-infinite plane 
approximation in the previous chapters. These are briefly 
discussed in the following. 

Shank et al suggested three design techniques in their 
paper [11] . 

1. The use of a series of stability-preserving transforma- 
tions to produce a stable two-dimensional recursive discrete 
filters from a stable one-dimensional continuous filter. 

2. A space domain design technique in which the unit 
sample response of the filter being designed is made to 
approximate a given desired unit sample response. 

3. A stabilization technique in which the denominator 
polynomial is replaced by its double planar least-squares 
inverse. Planar least-squares can be defined as follows: 

We have some given array C. We would like to find 
an array P such that C convolved with P approximates the 
unit pulse array U. That is, 

C * P « U 

122 






where the symbol * denotes two-dimensional convolution. In 
general it will not be possible to make C * P exactly equal 
to U. In actuality, C * P will be equal to some other 
array, G. If we choose P such that the sum of the squares 
of the elements of U-G is minimized, P is called a planar 
least square inverse of C. Planar least square inverse of 
P is called double planar least square inverse [11] . This 

new polynomial is conjuctured to have no zeros on 

-2 A , , , , 

u - {z^z^: | ZjJ ^_, |z 2 | >_ 1} and to have a frequency 

response which closely approximates that of the original 
transfer function. 

The first technique is known as a filter design with 
rotating axes and these filters are called rotated filters. 
This design procedure consists in considering a transfer 
function in the analog frequency domain, rotating the axes, 
and applying the bilinear transformation to obtain a two- 
dimensional recursive filter. In general, the starting 
point is a separable filter and stability of the resulting 
filter is not guaranteed even if the prototype filter is 
stable. This technique was used by Costa and Venetsanopoulus 
[23] to design low-pass filters having circular symmetry, 
and also they proved that rotated filters are stable if the 
angle of rotation is between 270° and 360°. A technique 
similar in character but rotated axes in digital frequency 
domain ((z,,z 2 ) domain) has been given by Chang and Aggarwal 
[24]. Hirona and Aggarwal [25] have accomplished the rota- 
tion through the introduction of rational powers of z^ and 



123 






z 2 , and the interpolation of input and output sequences. 
They have proved that, in this technique, the resulting 
filter is stable if the prototype filter is stable. Further, 
the distortion, a consequence of the bilinear transformation, 
is absent in the resulting filter when compared with Shank's 
et al. 's filter. 

The authors of [11] were unable to prove that the new 
denominator polynomial produced by the third technique, 
which is an extension of the one-dimensional case, was 
stable in the two-dimensional case. Nonetheless, several 
hundred filters have been designed using this technique 
without a counter example appearing [18]; furthermore, Jury, 
Kolavennu and Anderson [26] have shown it to be valid for 
certain low order cases. However, Genin and Kamp [27] have 
recently shown that counter-examples do, in fact, exist. 

Spectral factorization is another technique which can be 
used to produce a stable filter from an unstable one. By 
spectral factorization we mean a complex map that carries a 
stable rational transfer function into another stable transfer 
function exhibiting a different frequency response, at the 
same time maintaining some desirable characteristic [29] . 
Pendergrass, Mitra and Jury [28] used stability-preserving 
spectral transformations to modify the frequency response of 
quarter-plane filters. 

Mitra and Chakrabarti [29] extended this work. Read 
and Treitel [30] suggested a two-dimensional Hilbert trans- 
form to perform the spectral factorization, while Pister 



124 



[31] and Ekstrom and Woods [32] made use of the two- 
dimensional cepstrum. Ekstrom and Woods also suggested ways 
that the factors could be truncated and smoothed to preserve 
stability. They also showed that the usual quarter-plane 
filter is not the most general two-dimensional recursive 
filter. It is instead possible to have a recursive filter 
which in the limit has a unit sample response with support 
on an asymmetric half-plane, and also they developed 
stability conditions for this new class of filters. 

B. TWO-DIMENSIONAL DIGITAL FILTER DESIGN PROCEDURE BY 
DOUBLE Z-TRANSFORM 

1. Double Z-transform 

In the one-dimensional recursive digital filter 
design, the bilinear z-transform is often applied to an 
analog filter transfer function. We proposed that this one- 
dimensional design approach can be extended to the two- 
dimensional case. Before going into the design procedure, 
two-dimensional recursive filter design via the double 
bilinear z-transform, it is better to prove that we can 
use double z-transform for two-dimensional systems. 

Each of the canonic factors, proposed in Chapter 
IV, either separable or non-separable, can be converted 
into discrete equations by the double z-transform, as 
follows: 

2 1 - Z I' 



125 



2 X " z 2 
1 + z 2 



The single z-transform is equivalent to trapezoidal 
integration which can be shown as follows: 

X(s) = S Y(s) (5.2) 



.<*> = &J§i. 



or 

n 

Y n ~ y n _ x = / x(t) dt 
n-1 



?(x + x .) (5.4) 

2 n n-1 



where T is the sampling period. 

We write (5.4) in the z transform domain as 



(1 - z 1 )Y(z) = |d + z X )X(z) 



Hence 



X(z) = §( 1 " \[ ) Y(z) (5.5) 

1 + z 



126 



Comparing (5.5) with (5.2), it is seen that 

2 1 - z" 1 
S = J(i ^- T ) (5.6) 

1 + 2 

Equation (5.6) is known as a trapezoidal bilinear z -trans form, 

We can drive the double z-transform in a similar 
way, as follows: 



or 



Then 



,2 
W(x,y) = 33^ Z(x,y) (5.7a) 



W(s,,s 2 ) = s i s 2 z ( s ]_' s 2^ (5.7b) 



/ /d Z(x,y) = / /W(.x,y) dx, dy (5.8) 



by a discrete approximation, the left-hand side of (5.8) 
becomes (see Fig. 5.1) 



/ /d 2 Z(x,y) = [Z(x n ,y n " z ( x n' y n-l )] 

" ^V-l'V " Z(x n-l' Y n-l )] 
-1 -1 . -1 -1, 



(1 - z x - z 2 + z 1 z 2 )Z(.z 1 ,z 2 ) 



(5.9) 



127 



and the right-hand side of (5.8) is approximated as follows: 



//W(x,y)dxdy = A^£ [w(Xn ,y n ) + W(x n ,y n-1 ) 



+ "(Vl'V + W(x n-l'*n-l )] 



AXAVr, -1 , -1 , "I -1, TT/ V 

= — -^- [1 + z, + z 2 + z i z 2 ]W(z 1 ,z 2 ) 



(5.10) 



By equation (5.10) and (5.9) we have 



W(z 1 ,z 2 ) 1 - z" 1 - z" 1 + z~ z 2 



Z(z,,z ) . -1 -1 -1 -1 AxAy 
1 2 1 + z, + z 2 + z, z 2 2 



(5.11) 



or 



-1 i - 1 

2 1 " Z l 2 X ~ Z 2 

= ^ ( T—*^ ) ^ ( TT7^ ) (5 * 12) 

1 + z, J 1 + z 2 



By comparing (5.1a) and (5.2b) it is seen that 

2 1 " z l 

1 " Z 2 

Q.E.D. 

Thus, this is called a double z-transform. The 
double z-transform seems to be a reasonable substitution 



128 



y(t) 




x (O 




— n-1 n 



► t 



(a) One-dimensional case 



Z(*,y) 




w(x,y) 



>W(x,y) 




(b) Two-dimensional case 



Fig. 5-1- Two-dimensional extension of the 
Bilinear Transformation 



129 



for each canonic factor. However, each factor must be 
checked for stability, which is discussed in detail later in 
this chapter. If each transformed factor is BIBO stable the 
overall transfer function will be BIBO stable. But sometimes the 
double z-transform may in certain cases lead to unstable 
solutions [33] . Finally, it is noted that the double z-transform 
results in a quarter-plane symmetry which is a severe limitation. 
2 . Design Procedure 

The two-dimensional recursive digital filter design 
investigated in this thesis is the double bilinear z-transform. 
Proceed as follows: 

1. Determine a two-dimensional analog transfer 
function, 

A(s 1 ,s 9 ) 

hi«i'-2> = niTTifr (5 - 12> 

where A(s,,s 2 ) and B(s,,s 2 ) are mutually prime two-dimensional 
polynomials and B(s 1 ,s 2 ) satisfies 

B(s 1 ,s 2 ) f 

(5.13) 

V(s 1 ,s 2 ) z r = {( Sl ,s 2 ): Re (s 1 ) >_0,Re(s 2 ) >_ } 

in such a way that, the frequency response of G(s 1 ,s 2 ) 
meets the given specifications (for example as discussed 
in Chapters III and IV) . 



130 



2. Apply the double bilinear z-transform 

2 1 " z l 1 

5 1 = A { ^T } (5.14a) 

1 + z, 

2 1 " z l 

5 2 = -(- ij) (5.14b) 

1 + z 2 

where A is the sampling period, which is assumed the same in both 
directions. This yields the two-dimensional digital transfer 
function 



2 1 -z, 1-2, C(z,,z ) 

H( f ( TI ) 'f ( =T )} = T ^V Z 2 ] ~ D(z z ) (5 ' 15) 

A 1 + z 2 X A 1 + Z2 1 X 2 D( Z;l ,z 2 j 



where C(z,,z 2 ) and D(z,,z 2 ) are mutually prime two-dimensional 
polynomials in z, and z 2 . (5.15) will have the desired 
frequency response with proper sampling time without 
aliasing effect. Since 



, 1 - z? 1 

si = !< *r 

1 A 1 + z. 



maps the closed unit disc onto the closed right half 
plane, one would expect (5.15) to satisfy the following 
sufficient stability conditions [11], [20]. 



D(z 1 ,z 2 ) ? 



131 



\AW £ u2 - {(.z v z 2 ): | Zl | >_ 1, 



z 2 | >_ 1} (5.16) 



Unfortunately, however, the conditions (5.13) and (5.16) 
as we will show later in this chapter, are not always 
satisfied. 

In summary, we have two conditions to be satisfied 
for stability. 



I. B(s lf s 2 ) ? o 



\/(s 1 ,s 2 ) E r = {(s 1 ,s 2 ): R e s 1 >_ 0,R e s 2 > 0} (5.17) 



II. D(2 1 ,z 2 ) ? 



Y(z 1 ,z 2 ) £ u = {(z 1 ,z 2 ): Iz-J > l/|z 2 | >_ 1} 

(5.18) 
where V stands for "for each", and u is the closed unit bidisc 

C. DIGITAL CANONIC FACTORS AND STABILITY CONSIDERATIONS 
1 . Introduction 

Certain properties of two variable polynomials and 
rational functions which will be needed later will now 
be discussed. It is well known that a two variable polynomial 
is not, in general, factorable into first-order polynomials; 
but, a two variable polynomial can be factored into irreduci- 
ble factors which are themselves two variable polynomials 
but which cannot be further factored. Of course a given 



132 



polynomial may itself be irreducible. These irreducible 
polynomials are unique up to multiplicative constants [33] . 
Two polynomials which have no irreducible factors in common 
are said to be mutually prime. Consider two variable 
rational functions, in general, 



P (z, ,z 9 
G(z, ,z ) = 



l'"2' Q(z 1 ,z 2 ) 



where P(z^,z 2 ) and Q(z,,z 2 ) are mutually prime. Using 
terminology of [34], a point (zj,zl) making Q(z|,z 2 ) = 
but P(z',z') ^ is called a pole or a nonessential singu- 
larity of the first kind (such a point is analogous to a 
pole in the one variable case). A point (z£,z£) making 
P(zV,z~) = Q(zV,z 2 ) = is called a nonessential singu- 
larity of the second kind (such a point has no analogue in 
the one variable case). Clearly, if (zi,zl) is a pole, 
G(z,',z') = °°, but if (zV,z") is a nonessential singularity 
of the second kind, then G(zV,z") is undefined. 
2. Separable Factors 

The most common canonic separable factor in Chapter 
IV is 



H(<jj,,u) ) — "7t — —j, \ 1 1 . a r (5.19a) 

1 2 (1 + Joj, x,) (1 + 3oj 2 t 2 ) 



assuming x = x_ = 1. Substituting 



jw-L - S 1 , 



133 



ju> 2 = S 2 



We get 



H(s i' s 2» " (TTT[flTTip < 5 - 19b > 

Substituting equations (5.1a and 5.1b) into (5.19b), 
performs the double z-transformation. Thus, 

8(1 + z" 1 ) (1 + z" 1 ) 

T(z r z 2 ) = i- £_ (5.20) 

(a + z 1 ) (a + z 2 ) 



where 



A 2 



(A -2) 2 



a = 



a + 2 
a - 2 



A = sampling period 

This separable canonic form obviously satisfies the condition 
of (5.17). We should be careful about the second condition 
of (5.18) . 

When z, = z 2 = -1, we have a nonessential singularity 
of the second kind if a = 1, but a = 1 never occurs. 
Therefore, there is no nonessential singularity of the second 
kind. Thus, this separable case is stable by satisfying 



134 



two essential necessities (5.17) and (5.18). The frequency 
response of this separable case is plotted in Fig. 5.2, with 
sampling time A = 0.08 which fits that of the analog frequency 
response. 

The double z-transform for equation (5.19a) is 



A 2 (l + z 1 1 ) (l + z^ 1 ) 



T(Z 1 /Z 2 ) = Ta-2 T;l ) (A-2t 2 ) A+2 T;l ~ aT2~^ ~ 

( A^2T7 +Z 1 )( A^2TJ +Z 2 > 

(5.21) 

and this general case for the two-dimensional separable 

factor is stable by satisfying the stability necessities 

of (5.17) and (5.18) 

The second analog canonic factor in Chapter IV is 



H(j Ul rJ0. 2 ) (1 + jo) 1 T 1 ) (1 + JW 2 T 2 ) (1 + j w 2 T 3 ) ' 

When we perform the double bilinear transform for 
this, we get 



3 -1 -12 
A J (1 + z/) (1 + z 9 x r 

T(z n ,z„) = 



1' 2' A+2t, , A+2T- , A+2t-. , 

(A-2x 1 )(A-2x 2 )(A-2x 3 )[^ + z- i ][ I=2 ^ + z 2 H I=n i+«; ] 

(5.22) 

The frequency response of (5.22) for t, = 1, t 2 = 2, t^ = 4 
is shown in Fig. 5.3 which almost fits the analog case. 



135 





Fig. 5.2a. The frequency response of 



T(z 1 ,z 2 ) = 



where 3 = 



e(l+z 1 " 1 ) (l+z 2 1 ) 

(a+z^ ) (a+z~ ) 

A 2 A+2 
~ , a = 



(A-2)" 

(dB vs. log a), , log u}~) 



A-2 



= 0.08 



136 




Fig. 5.2b. The frequency response of 



T(z 1 ,z 2 ) = 



6(l+z 1 ~ 1 ) (l+z 2 1 ) 
(a+z 1 " 1 ) (a+z 2 -1 ) 



3 = 



A' 



(A-2) 



2 ' 



a = |i| , A = 0.08 



(dB vs. a). ,0)-) (linear frequency scale) 



137 



IT). 



0. -, 



*4 



10.-, 




Fig. 5.2c. 



The contour plot of the frequency response of 
3(l+z 1 " 1 ) (l+z 2 _1 ) 



T(z x ,z 2 ) 



6 = 



2 ' 



(a+z 1 ~ ) (a+z 2 ) 
a = f±§, A= O.OB 



where 



(A-2) 
(linear frequency scale) 



138 



We can summarize the m xn analog separable factor as 



M N 
H(ju> ,ju, ) = [ IT , .-*• ] [ n , , / ] (5.23) 



After performing double z-transform, (5.23) becomes 



M -1 N -1 

M a (1 + z, ) N A (1 + z,) 

T(z,,z ) = [ n — =i ][ n — | ] 

(A-2t.) (-— ii+zT 1 ) 3 L (A-2x.) (_-^J-+z: 1 ) 
1 A-2x -1 3 A-2t ■ 2. 

15.24) 

Equations (5.23) and (5.24) satisfy the conditions 
of (5.17) and (5.18). As a result, a separable filter 
results in separable filters in digital domain and they are 
stable. 

3. Non-separable Factors 

The first non-separable factor in Chapter IV is 
4.d with a = -1, p = q = t, = t 2 = 1 as given below 



1 



H ( joj, t j oj ) - j— - ■ 

- 1 1 2 1 + 3oj-i](jo~ 



(5.25a) 



or 



«<V S 2> ■ TTT^ (5 - 25b: 



139 




Fig. 5.3a. The frequency response 

A 3 (l+z 1 " 1 ) d+z^ 1 ) 
T(z 1/ z 2 ) = 



(A-2) (A-4) (A-S, (f±| .z,- 1 ) (f±| .zf 1 ) <£±f +z 2 -h 



where A = 0.0 8 
(dB vs. log co, , log u^) 



140 




Fig. 5.3b. The frequency response of 



T( 2l ,z 2 ) = 



A 3 (l+z 1 " 1 ) U+z^ 1 ) 



(A-2) (A-4) (A-8) (|±| +z 1 " 1 ) (|±i +z x X ) (f±f + z 2 1 ) 



where A = 0.08 



(dB vs. u.,ii)J 



141 




!S.u 



Fig. 5.3c. The contour plot of the frequency response of 



T(z 1 ,z 2 ) = 



A 3 (l + z 1 " 1 ) (l+z^ 1 ) 



(A-2) (A-4) (A-8) (|±| +z 1 - 1 ) (|±J +a 1 " 1 ) <£±§ -hz^ 1 ) 



where A = 0.0 8 
(linear frequency scale) 



142 



and by the double 2 -trans form, we get 



A 2 (l +Z, 1 ) (1 + Z' 1 ) 

T(Z 1' Z 2 ) = — 2 =T TT T Tt-ZT (5 ' 26) 

(A +4) (A -4) (z L + z l ) + (A +4)z 1 z X 



In the analog domain (5.25b) does not satisfy the 
condition (5.17). We have a singularity at 



Ul u) 2 = 1 



For the second condition, we have a nonessential 
singularity of the second kind at z, = -1, z_ = -1. This 
does not tell us whether it is stable or not, When we 
check the frequency response, it is seen that this factor 
is unstable. 

The second non-separable canonic factor can be 
derived from (4.d) with a = q = -1 and T-,=T 2 = las 
follows. 



H(ju)i ,jui ) = 3 (5.27a) 

1 + 



2 tu. 



or 



H < s i< s 2> - — 4r = sr^r < 5 - 27b > 

1 + IT 
s 2 



143 



Then, (5.27b) becomes after performing the double z-transform, 
as follows: 

(l + zT 1 ) (1 -z' 1 ) 

T(z ,z ) = _ 2 ( 5.28) 

2U-« 1 r « 2 I ) 



This canonic factor does not satisfy either 
stability condition. Hence it is unstable. 

The other series of canonic factors can be derived 
from (4.e) . 

For example, with a = -1, t 1 = t- « p » 1, we 

get 



HtjoiwjwJ = -, , - 1 , ■ (5.29a; 

-L ^ 1 + 3 co, +3 u^ 



or 



H(S 1' S 2» = l + s^ + s 2 ( 5 - 29b> 



or in the digital domain 



A(l + z" 1 ) (1 + z" 1 ) 



T(z,,z ) = — —, _ (5.30) 

(A+4) +A(z x x + z 2 ) + (A-4)z 1 z 2 



This canonic factor satisfies the first condition but 
it has a nonessential singularity of the second kind, at 
z, = z~ = -1 and so we should check the frequency 



144 



response. We see that there is a singularity at 

(u) 1 ,o) 2 ) = (it, it). Thus this non-separable canonic factor 

is unstable. 

4 . Conclusions 

All of the canonic factors discussed in this 
chapter are collected in Table I. By comparing the results 
shown in this table, we can say that separable stable canonic 
factors result in stable separable filters. As seen from 
Table I, the separable transfer functions 1, 2 and 3 which 
start with stable separable ones in the two-dimensional 
continuous frequency domain end up as stable separable trans- 
fer functions in the two-dimensional discrete domain. This 
is also true for the mxn separable canonic factor 4. 

Non-separable canonic factors have very different 
characteristics. The unstable non-separable case ends up 
with a non-separable case in the discrete domain, which is 
also unstable. Although we could not prove this result 
generally, we also could not find any counter example to 
exist. The transfer functions of 5, and 6 start in the 
unstable case, and end up with the unstable case in the 
discrete domain. The transfer functions of 7 starts with 
the stable case in the continuous frequency domain, ends up 
with the unstable case in the discrete domain. Thus we 
have to be very careful with the double z-transform since it 
does not yield stable results every time. It may yield 
unstable results even if it starts from the stable case. 



145 






UJ 

_J 

CO 

< 



>- 



CO 
< 

co 



C\J 
M 



CM 
<f) 

J) 



LU 



GQ 
< 









ea 












H 









LU 



CD 
< 



CO 






I m 



w 






I 



NT 
I 



4- 






CM 



LU 



co 
< 



uo 



J-flL 



If 









<< 



V 



I 

i 
X 

I 









00 

+ 



rO 



LU 



CD 
< 



CO 






r 


r ^-^ \ 




V 




♦ 


*" s 




t-r 






1 

< 


~ • 


*~<+ 


«- ^ 


.**• 




C«i 


*< 








II 



rJ 

WO 



*=* 



(/T 

+ 



i ! i 



«■ 



3iaVd3d 3S 



146 



LU 

CD 
< 






LU 



CQ 
< 



UJ 

—I 

cn 

< 



ui 



i « 

4- 



i — 

m 

+• 



-< 



fW 



J" 









«i 






1 H 




(N 


*— «. 


\ 


1 m 


— 


hi 


-.jj 


• 




1 ^ 


^ 


r« 


i . 


i 


fV 


— 


+ 




_ 


N 


'-" / 





I .4 






> * 

"7. 



» 



• -. 
I* 

ru 

+ 



+ 



*W 



iT) 



CD 



+ 

v/T 

4» 



1^ 



319VyVd3SNn 



147 



VI. CONCLUSIONS 

A technique for the use of semi-infinite planes to 
approximate the log magnitude versus log frequency charac- 
teristics of two-dimensional filters is presented. This 
technique is applied to quarter plane filters and works 
well for separable transfer functions. Other non- 
separable canonic (basic) transfer functions are also 
studied in terms of their planar approximation. 

The following separable and nonseparable canonic factors 
have been studied. 

(1) (1 + jo) 1 x 1 ) p (l + jw 2 T 2 ) q 

(2) (1 + oo 1 x 1 ) P (l + jco 2 T 2 ) q (l + joj 2 T 3 ) q ' 

(3) [1 + (jco 1 t 1 ) P ] a 



q,a 



(4) [1 + (joo 1 T 1 )^(l + ja> 2 T 2 )*] 



p,a 



(5) [1 + (ja) 1 T 1 + jco 2 T 2 )^] 



The properties and stability of these canonic factors 
are investigated first in the continuous variable domain 
and then after the double bilinear z-transform is applied 
to extend them to two-dimensional recursive digital filters 
These are all summarized in Table I in Chapter V. An 
important result is that although the double bilinear 
z-transform works well for separable factors, it does not 



148 



always lead to a stable digital realization starting with a 
stable analog transfer function. 

The following are suggested for further research: 

1. Other types of canonic factors need to be considered. 
For example, the general form: 

H(j» lf jw 2 ) = [1 + (ja3 1 T 1 ) P +(ja J2 T 2 ) q +(jco 1 T 3 ) r (jw 2 T 4 ) t ] a 

2. The work should be extended to non-causal (half- 
plane) filters. 

3. Transfer functions which have circular symmetry in 
the discrete domain need to be investigated. 

4. Another approach is to consider half-plane 
recursive algorithms, by working backwards from the discrete 
domain (z, ,zj to continuous domain (s,,s 2 ). 

5. Stability of canonic factors needs to be investigated 
in detail . 

6. The use of predistortion should also be considered. 



149 



LIST OF REFERENCES 



1. Oppenheim, A.V. and Schafer, R.W., Digital Signal 
Processing , Prentice Hall, Englewood Cliffs, New Jersey, 
1975. 

2. Rabiner, L. and Gold, B. , Theory and Application of 
Digital Signal Processing , Prentice Hall, Englewood 
Cliffs, New Jersey, 1975. 

3. Rabiner, L.R., et al., "Terminology in Digital Signal 
Processing," IEEE Trans. Audio Electroacoust . , vol. AU-20, 
p. 322-334, December 1972. 

4. Stockham, T.G., "High Speed Convolution and Correlation," 
AFIPS Conference Proceedings, 28, pp. 229-233, 1966. 

5. Hall, E.L., "A Comparison of Computations for Spatial 
Frequency Filtering," Proc. IEEE (Special Issue on Digital 
Picture Processing), Vol. 60, pp. 887-891, July 1972. 

6. Huang, T.S., "Two-Dimensional Windows," IEEE Trans. Audio 
Electroacoust. (Corresp.), Vol. AU-20, pp. 88-89, March 
1972. 

7. Hu, J.V. and Rabiner, L.R., "Design Techniques for Two- 
Dimensional Digital Filters," IEEE Trans. Audio Electro- 
acoust. (Special Issue on Digital Filtering), Vol. AU-20, 
pp. 249-257, October 1972. 

8. Arazi, E., "Two-Dimensional Digital Processing of One- 
Dimensional Signal," IEEE Trans. Acoustics, Speech, and 
Signal Processing, Vol. ASSP-22, pp. 81-86, April 1974. 

9. Kamp, Y. and Thiran, J. P., "Maximally Flat Nonrecursive 
Two-dimensional Digital Filters," IEEE Trans. Circuits and 
Systems , Vol. CAS-21, pp. 437-449, May 1974. 

10. Rabiner, L.R., McClellan, J.H., and Parks, T.W. , "FIR 
Digital Filter Design Techniques Using Weighted Chebyshev 
Approximation," Proc. IEEE , Vol. 63, pp. 595-609, April 
1975. 

11. Shanks, J.L., Treitel, S. and Justice, J.H., "Stability 
and Synthesis of Two-Dimensional Recursive Filters," 
IEEE Trans. Audio Electroacoust. , Vol AU-20, pp. 115-128, 
June 1972. 



150 



12. Mitra, S.K., Sagar, A.D., and Pendergrass, N.A., 
"Realizations of Two- Dimensional Recursive Digital 
Filters," IEEE Trans. Circuits and Systems , Vol. CAS-22, 
pp. 177-184, March 1975. 

13. Green, B.W., and O'Handley, A.D., "Recent Developments 
in Digital Image Processing at the Image Processing 
Laboratory at the Jet Propulsion Laboratory," Proc. IEEE , 
Vol. 60, pp. 821-827, July 1972. 

14. Treitel, S., Shanks, J.L. , and Frasier, C.W., "Some 
Aspects of Far Filtering," Geophysics , Vol. 32, 

pp. 789-800, October 1967. 

15. Huang, T.S., Schreiber, F.W., and Tretiak, O.J., "Image 
Processing," Proc. IEEE , Vol. 59, pp. 1586-1609, 
November 1971. 

16. Hunt, B.R., "Digital Image Processing," Proc. IEEE , 
Vol. 63, pp. 693-708, April 1975. 

17. Wood, L.C. and Treitel, S., "Seismic Signal Processing," 
Proc. IEEE , Vol. 63, pp. 649-661, April 1975. 

18. Mersereau, Russell M. and Dudgeon, D.E., "Two- 
Dimensional Digital Filtering," Proc. IEEE , Vol. 63, 
pp. 610-623, April 1975. 

19. Chong, H. and Aggarwal , J.K., "Fundamentals of Two- 
Dimensional Recursive Filtering," Digital Signal 
Processing , (Edited by Aggarwal, J.K.), pp. 211-260, 
Pt. Lobos Press, CA. , 1979. 

20. Huang, T.S., "Stability of Two-Dimensional Recursive 
Filters," IEEE Trans. Audio Electroacoust . , Vol. AU-20, 
pp. 158-163, June 1972. 

21. Ansell, H.G., "On Certain Two-Variable Generalizations 
of Circuit Theory, with Applications to Network of 
Transmission Lines and Lumped Reactances," IEEE Trans. 
Circuit Theory , VoL CT-11, pp. 214-223, June 1964. 

22. Anderson, B.D.O. and Jury, E.I., "Stability Test for 
Two-Dimensional Recursive Filters," IEEE Trans. Audio 
Electroacoust . , Vol. AU-21, pp. 366-372, August 1973. 

23. Costa, J.M., and Venetsanopoulos , A., "Design of 
Circularly Symmetric Two-Dimensional Recursive Filters, 1 
IEEE Trans. Acoustics, Speech and Signal Processing , 
Vol. ASSP-22, pp. 432-443, December 1974. 



151 



24. Chang, H. and Aggarwal, J.K., "Design of Two-Dimensional 
Recursive Filters by Interpolation/' IEEE Trans. Circuits 
and Systems , Vol. CAS-24, pp. 281-291, June 1977. 

25. Hirono, K. and Aggarwal, J.K., "Design of Two- 
Dimensional Recursive Digital Filters," IEEE Trans. 
Circuits and Systems , Vol. CAS-25, pp. 1066-1076, 
December 1978. 

26. Jury, I., Kolavennu, V.R., and Anderson, B.D.O., 
"Stabilization of Certain Two-Dimensional Recursive 
Digital Filters," Proc. IEEE , Vol. 65, pp. 887-891, 
June 1977. 

27. Genin, Y.V. and Kamp , Y.G., "Two-Dimensional Stability 
and Orthogonal Polynomials on the Hypercircle, " Proc. 
IEEE , Vol. 65, pp. 873-881, June 1977. 

28. Pendergrass, N.A., Mitra, S.K., and Jury, E.I., "Spectral 
Transformations for Two-Dimensional Digital Filters", 
IEEE Trans. Circuits and Systems , Vol. CAS-23, pp. 26- 
35, January 197 6. 

29. Chakrabarti, S. and Mitra, S.K., "Design of Two-Dimensional 
Digital Filters via Spectral Transformations", Proc. 

IEEE , Vol. 65, pp. 905-914, June 1977. 

30. Read, Randol. R. and Treitel , S., "The Stabilization of 
Two-Dimensional Recursive Filters via the Discrete 
Hilbert Transform" , IEEE Trans. Geoscience Elect . 

Vol. GE-11, pp. 153-160. July 1973. 

31. Pister, P., "Stability Criterion for Recursive Filters", 
IBM Journal of Research and Development , Vol. 18, 

pp. 59-71, January 1974. 

32. Ekstrom, P.M. and Woods, J.W., "Two-Dimensional Spectral 
Factorization with Applications in Recursive Digital 
Filtering", IEEE Trans. Acoustics, Speech and Signal 
Processing , Vol. ASSP-24, Pp. 115-128, April 1976. 

33. Goodman, D., "Some Difficulties with the Double 
Bilinear Transformation in Two-Dimensional Recursive 
Filter Design", Proc. IEEE , Vol. 66, pp. 796-797, 
July 1978. 

34. Goodman, D., "Some Stability Properties of Two-Dimensional 
Linear Swift-Invariant Digital Filters", IEEE Trans . 
Circuits and Systems , Vol. CAS-24, pp. 201-208, April 
1977. 



152 



INITIAL DISTRIBUTION LIST 



No. Copies 



1. Defense Documentation Center 2 
Cameron Station 

Alexandria, Virginia 22314 

2. Library, Code 014 2 2 
Naval Postgraduate School 

Monterey, California 93940 

3. Department Chairman, Code 62 1 
Department of Electrical Engineering 

Naval Postgraduate School 
Monterey, California 93940 

4. Professor S. R. Parker, Code 62Px 10 
Department of Electrical Engineering 

Naval Postgraduate School 
Monterey, California 93940 

5. Deniz Kuvvetleri K-ligi 1 
Personel Sube Bsk-ligi 

ANKARA/TURKEY 

6. Deniz Harb Okulu K-ligi 1 
Heybeliada 

ISTANBUL/TURKEY 

7. Coskun Zumrutkaya 3 
Yuzbasiler cad. No=20 

Degirmendere 
KOCAELI/TURKEY 



153 



I<2 AU^Qo 



25787. 



Thesis 

Z8 Zumrutkaya ] 8381 6 

Geometric design tech- 
nique for two-dimensional 
analog and digital filters. 



c.l 



l,c A^jdO 



- 25 7 8 7^ 



Thesis 

Z8 Zumrutkaya 

c.l Geometric design tech- 

nique for two-dimensional 
analog and digital filters 



183816 



thesZ8 

Geometric design technique for two-dimen 




3 2768 001 07029 5 

DUDLEY KNOX LIBRARY